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ABSTRACT 

Angular momentum transport and accretion in protoplanetary discs are generally believed to 
be driven by MHD turbulence via the magneto-rotational instability (MRI). The dynamics of 
solid bodies embedded in such discs (dust grains, boulders, planetesimals and planets) may 
be strongly affected by the turbulence, such that the formation pathways for planetary systems 
are determined in part by the strength and spatial distribution of the turbulent flow. 

We examine the dynamics of planetesimals, with radii between 1 m - 10 km, embedded in 
turbulent protoplanetary discs, using three dimensional MHD simulations. The planetesimals 
experience gas drag and stochastic gravitational forces due to the turbulent disc. We use, and 
compare the results from, local shearing box simulations and global models in this study. 

The main aims of this work are to examine: the growth, and possible saturation, of the 
velocity dispersion of embedded planetesimals as a function of their size and disc parameters; 
the rate of radial migration and diffusion of planetesimals; the conditions under which the 
results from shearing box and global simulations agree. 

We find good agreement between local and global simulations when shearing boxes of 
dimension 4H x 16H x 2H are used (H being the local scale height). The magnitude of the 
density fluctuations obtained is sensitive to the box size, due to the excitation and propagation 
of spiral density waves. This affects the stochastic forcing experienced by planetesimals. The 
correlation time associated with the stochastic forcing is also found to be a function of the box 
size and aspect ratio. 

The equilibrium radial velocity dispersion, cr(y^), obtained depends on the radii, Rp, of 
the planetesimals. Bodies with Rj, = 50m achieve the smallest value with cr{Vf) ^ 20ms"'. 
Smaller bodies are tightly coupled to the gas, and boulders with /?p = 1 m attain a value of 
o"(vr) similar to the turbulent velocity of the gas (~ 100 ms '). Equilibrium values of o-(vr) 
for bodies larger than 100 m are not achieved in our simulations, but in all models we find 
rapid growth of the velocity dispersion for planetesimals of size 1 km and 10 km, such that 
cr(vv) > 160ms"' after a run time of 1200 orbits at a distance of 5 au from the central star. 
These values are too large to allow for the runaway growth of planetesimals, and mutual 
collisions would lead to catastrophic disruption. Radial migration due to gas drag is observed 
for bodies with R^ ^ Im, and is only modestly affected by the turbulence. Larger bodies 
undergo a random walk in their semi-major axes, leading to radial diffusion through the disc. 
For our fiducial disc model, we estimate that radial diffusion across a distance of ^ 2.5 au 
would occur for typical planetesimals in a swarm located at 5 au over a disc life-time of 5 Myr. 
Radial diffusion of this magnitude appears to be inconsistent with Solar System constraints. 

Our models show that fully developed MHD turbulence in protoplanetary discs would 
have a destructive effect on embedded planetesimals. Relatively low levels of turbulence are 
required for traditional models of planetesimal accretion to operate, this being consistent with 
the existence of a dead zone in protoplanetary discs. 

Key words: accretion disks - magnetohydrodynamics (MHD) - methods: numerical - plan- 
etary systems: formation - planetary systems: protoplanetary disks 
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1 INTRODUCTION 

The dynamical and collisional evolution of planetesimals is a fun- 
damental issue which needs to be understood if progress is to 
be made in developing a theory of planetary system formation. 
According to the core accretion theory, a process which begins 
with the collision and sticking of small dust grains within a pro- 
toplanetary disc lea ds eventually to the forma tion of kilometre 
sized planetesimals JWetherill & StewartI Il993h . Although alter- 
nativ e scenarios have been put forward for pl anetesimal forma- 
tion JGoldreich & War j[T973l : Ijohansen et"aij |2007). which avoid 
the requirement for such large bodies to grow via simple two- 
body agglomeration processes, the further growth of planetesimals 
into planetary embryos and cores generally requires planetesimals 
themselves to accrete via mutual collisions through a process of 
runaway growth, followed by oligarchic growth ( Ida & Makino 
ri99ilKokubo & IdJl998l) . 

Rapid runaway growth requires that the velocity dispersion 
of the planetesimal swarm remains significantly smaller than the 
escape velocity from the surfaces of the largest accreting planetes- 
imals, ensuring that gravitational focusing is important. For bod- 
ies of radius 10 km, and with internal densities Qp == 2gcm"^, 
the escape velocity is 10 ms"', and scales linearly with radius. 
Clearly this is a stringent requirement, which is easily met within 
a self-stirring planetesimal swarm whose size distribution is rea- 
sonably shallow, but which may be difficult to satisfy in the pres- 
ence of an external source of stirring. One such source may be 
turbulence within the protoplanetary disc. Planetary growth times 
which rely on mutual collisions between planetesimals occurring 
at rates which are determined by the geometric cross section are 
prohibitively long, leading to estimated planetary growth times 
which are much in excess of typical protostellar disc life times 
JHaisch et aOlOOlh . 

A further constraint during the runaway growth phase is that 
collisional velocities should be small enough to avoid catastrophic 
disruption of planetesimals. For bodies in the 1 - 10 km size range, 
for which self-gravity starts to become more important than mate- 
rial strength in holding planetesimals together, collisions between 
similar sized bodies with impact speeds which are modestly in ex- 
cess of the escape velocity will lead to bre ak up of the planetesimals 
rather that accretion and growth. Indeed Benz & Asphauj l l 19991) 
suggest that mutual collisions between 1 km sized bodies will result 
in catastrophic disruption if the impact speeds exceed ~ 20ms"', 
depending o n the material compos i tion of the impactors. In a more 
recent study, I Stewart & Leinhardtl ( 120091) suggest reduced impact 
speeds of ~ 10 ms"' will be destructive. Once again, we see that 
an external source of planetesimal stirring may prevent the rapid 
growth of planetary mass bodies by the accretion of planetesimals. 

The canonical mass accretion rate for T Tauri stars is ~ 
10"** Mq yr"' ( Sicilia-Aguilar et al. 2004). Such accretion rates re- 
quire a source of anomalous disc viscosity and angular momen- 
tum transport to operate, generally thought to be turbulence. The 
most likely sour ce of disc turbulence is the magneto-rotational 
instability (MRI. Is^bus & Hawlevlll991h which has been shown 
to develop into full non linear MHD turbulence in numerous 
studie s, using both local shearin g box simulation s (Hawlev et al. 
19951). and global simulat ions ( lArmitagel 1 19981 : fHawlev. ,2001. ; 



Papaloizou & Nelson 120031) . The nature and saturation state of 



MHD turbulence gen erated by the MRI is the subject of on-going 
study JFromang & P apaloizou 2007; Fromang et al. 2007). In this 
work, we use the dependence of the turbulent stresses and density 
fluctuation amplitude on the strength of the net component of the 



magnetic field to examine the evolution of planetesimals in discs 
with different levels of turbulence. We use simple disc models, 
which neglect non-ideal MHD effects and vertical stratification. As 
such, this is the first in a series of papers in which we examine 
how turbulence affects the dynamics of planetesimals embedded in 
turbulent discs. Future papers will explore the effects of vertical 
stratification and dead zones. 

There have been numerou s studies of planets embedd ed in 
turbulent pro toplanetary discs. Nelson & PapaloizoiJ ( 120031) and 
[winters et alj 112003) examined the formation of gaps by lovian 
mass p lanets, and the migration torq ues exerted by the disc on the 
planet. iNelson & Papaloizoii ( 120041) performed global simulations 
of low mass planets embedded in turbulent discs. They showed 
that such bodies are subject to fluctuating torques which should 
induce stochastic migration, and suggested that this might provide 
a means of mitigating against the rapid type I migration expected 
to occur for low mass planets jWarall997h . iLaughlin et al.l J2004h 
published a similar study using analytical fits to MHD simulations, 
and reached similar conclusions. Papaloizou et al. (2004) presented 
results from both global and local shearing box simulations con- 
taining both high and low mass planets, and showed good agree- 
ment between the simulation set-ups for predicting the transition 
between linear and non linear disc respon se to th e presence of 
an embedded planet. In a follow-up paper, iNelsonI (2005) exam- 
ined the orbital evolution of low mass embedded planets, showing 
that over simulations run times of ^ 100 orbits, turbulence induces 
stochastic migration for planets in the range 1-10 Me, and induces 
the growth of orbital eccentricity. In more recent work. lOishi et al.l 
(2007) examined the stochastic forces experienced by planets in 
stratified disc model s with and withou t dead zones using shearing 
box simulations, and l Yang et alj J2009I) examined the orbital evolu- 
tion of swarms of test particles embedded in non stratified turbulent 
discs. 

A significant volume of related work has examined the influ- 
ence of disc turbulence on embedded planets and planetesimals us- 
ing prescriptions or simple models for the effects of turbulence, 
[lohnson et al. ( 2006) developed a Fokker-Planck description for the 
stochastic evolution of planets, and examined the survival proba- 
bilities of distributions of planets subject to type I migration and 
superposed stochast ic migration. A simila r study has been pub- 
lished recently by Adams &Blocd J2009l) . lOgihara e t aT ('2007|) 
have used N-body simulations plus a prescription for stochastic 
forcing to examine t he eflF ects of turbulence on terrestrial planet 
formation. Ilda et al.l <2008l) used a similar prescription of turbu- 
lent forcing and examined the growth of eccentricity for planetesi- 
mals, exploring in particular the possibility of reaching catastrophic 
disruption velocities. Adams et al. ( 2008) and Rein & Papaloizoul 
(2009) examined the stability of mean mot ion resonances for pairs 
of planets embedded in turbulent discs, and lBaruteau & Linl <201CI) 
have examined the saturation of corotation torques in turbulent 
discs by means of hydrodynamic simulations subject to turbulent 
stirring. 

In this paper, we examine in detail the orbital evolution of 
planetesimals of different size (ranging between 1 m and 10 km) 
embedded in turbulent disc models by means of three dimensional 
MHD simulations. A key issue that we explore is the set of con- 
ditions and numerical parameters that provide good agreement be- 
tween local shearing box simulations and global simulations. We 
find that it is possible to obtain good agreement between these two 
numerical set-ups, provided that the shearing box dimensions are 
chosen appropriately. Other important issues that we examine in- 
clude the growth of the velocity dispersion of embedded swarm of 
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planetesimals, and the saturation value of this velocity dispersion as 
a function of planetesimal size due to a balance being achieved be- 
tween gas drag and turbulent forcing. We explore the implications 
of our results for the efficiency of runaway growth of planetary em- 
bryos, and the possibility that planetesimals may enter a phase of 
catastrophic disruption through mutual collisions. We also exam- 
ine the rate at whic h planetesimals migrate due to both gas drag 
induced radial drift jWeidenschillinal 19771) . and due to diffusion 
caused by the fluctuating gravitational field of the turbulent disc. 
We examine under which conditions each of these processes are 
dominant, and we explore the implications of our results for the ra- 
dial drift of planetesimals in the Solar nebula and limits that might 
be placed on the magnitude of turbulent fluctuations which were 
present during the early phases of Solar System formation. 

This paper is organised as follows. In Sect|2]we describe the 
numerical set-up and the parameters of the disc models. In Sect.|3] 
we present our results. In Sect.|4]we discuss our work in the context 
of previous work and draw our conclusions. 



2 MODEL DESCRIPTION 

We perform self-consistent simulations of hydromagnetic turbu- 
lence using two diff^erent set-ups: shearing box simulations which 
represent a local patch of a protoplanetary disc; global disc models 
which simulate a larger section of a protoplanetary disc and include 
the full set of curvature terms in the equations of motion. A major 
goal of this work is to compare the results of these two different 
approaches. 

A key question that needs to be addressed is for what dimen- 
sions of the the shearing box, in units of the local scale height H, 
do density fluctuations created by the turbulence reach a converged 
amplitude and spectrum, and do these match the results of global 
models. To make the comparison as straightforward as possible, we 
neglect vertical stratification and assume an isothermal equation of 
state. 

In both configurations, the hydromagnetic turbulence is driven 
via the non-linear development of the magneto-rotational insta- 
bility. At present, the issue of the saturated amplitude of MRI 
turbulence remains unresolved and is a topic of ac tive research 
JFromang & Papaloizoj 120071 : iFromang et al.l 1200'^ . In the ab- 
sence of a better alternative, we therefore adopt a practical perspec- 
tive and impose a net vertical or azimuthal flux, for which numeri- 
cal convergence can be obtained (Davis et al. 2009). Negl ecting the 
depen dence on the magnetic Prandtl number l iLesur & Longarettil 
I2007I) . we furthermore restrict ourselves to the case of ideal MHD. 
This approach is justified by the observed correlation between the 
strength of the turbulence and the amplitude of the resulting den- 
sity fluctuations (Yangetal. 2009). This means that we regard 
the strength of the imposed field as a control parameter which 
can be tuned to vary the turbulence amplitude in the local and 
global context. The global cylindrical disc models are computed 
with a modified version of t he original finite difference code nir- 
vana jZiegler & Yorkelll997h . For the local shearing box models, 
we make u se of the new l y dev eloped second-order Godunov code 
NiRVANA-iii ( IZieglej2004l.l2008h . 



2.1 Numerical methods - local model 

For our standard model, we adopt a box siz^j of 4x 16x 2 pressure 
scale heights /f at a resolution of 32 grid points per H. Boundary 
conditions are periodic in the azimuthal (y) and vertical (z), and 
sheared-periodic in the radial (x) direction. The initial net vertical 
magnetic field corresponds to a plasma parameter j6 == 6000, result- 
ing in a typical saturation level a ^ 0.05 of the turbulence, wher e 
a is the effective viscosity parameter dShakura & Svunvaevll 19731) . 
andyS is ratio of thermal to magnetic pressure. 

Because the gas drag forces acting on massive particles de- 
pend on the actual physical value of the gas density, we have to 
prescribe a set of conversion factors to link our model to a repre- 
sentative protoplanetary disc. We chose a fiducial radius Rq = 5 au, 
and a geometric disc thickness of H/R = 0.05 at /? = 1 au. Note that 

this aspect ratio is scaled with i?''* to be consistent with the Hayashi 

I II 1 

minimum mass solar nebula (MMSN, Hayashi 1981), yielding a 

value of =^ 0.075 at 5 au. Furthermore, we chose a slightly higher 

average mass density than in this model to yield a column density 

S = 160 gcm"^ and sound speed Cs = 1 kms"' comparable to the 

global simulations. 

For the local model, we evolve the following set of non-linear 

partial differential equations 



a,(£.v) + v. 



B- BB 

gy\ + (p+—)l 

a,B-Vx(vxB) 



0, 

-2£)fl z X (v + qQ.xf) 
0, (1) 



comprising the standard formulation of ideal MHD in the shear- 
ing box approximation, and where we have assumed an isothermal 
equation of state p = gc^ and neglected the effects of stratification. 
The two momentum source terms are the Coriolis force -IClixy in 
the locally corotating frame, and the tidal term IqOrxx, with shear- 
parameter q = 3/2, describing the linearised effect of the Keplerian 
rotation|j Care has been taken in implementing t he source terms to 
minimise the error in the epicyclic mode energy ( Gressel & Zieglei| 
2007), albeit not to the extent w here it is conserved to machine ac- 
curacy JStone & GardinenuOlOn . 



2.1.1 Numerical scheme & orbital advection 



As has been recently demonstrated by ^B alsara & Meveg ( 120100 . 
the adequate modelling of the MRI with finite volume codes de- 
pends on the reconstruction strategy used and, in particular, on 
the ability of the Riemann solver to capture the Alfven mode. To 
improve the treatment of discontinuities in the Godunov scheme, 
we therefore extended nirvana-iii with the Harten-Lax-van Leer 
Discontinuities (H LLD) approximate Riemann solver proposed by 
iMivoshi & Kusanol ( |2005|) . 

In accordance with its finite volume approach, the nirvana-iii 
code evaluates the components of the electromotive force (EMF) 
at cell interfaces. Since the discretisation of the constrained trans- 
port algorithm intrinsically requires edge-aligned EMFs, some sort 
of interpolation is required. In its origin al form, nirvana-iii imple- 
ments the arithmetic avera ge proposed by Balsara & SDicen ( ll999h . 
As discussed in Sect. 3.2 of lGardiner & Stond ( 120051) . this approach 
however lacks the required directional biasing to guarantee the sta- 
bility of the numerical scheme. It has further been demonstrated by 



See Sect. l3.3l for a discussion on the effect of the box size. 
^ The variable x = R- Rqis the radial displacement from the box centre. 
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iFlock et alj 1 120091) that this can lead to artificial growth of instabil- 
ities in the context of net-flux MRI, and we have reproduced this 
result. To resolve this issue, we have successfull y implemented and 
tested the upwind reconstruction procedure of iGardiner & Stond 
( l2005h . 

Following the long-term evolution of a shearing flow in boxes 
of substantial radial extent puts high demands on computational re- 
sources. For a Keplerian rotation profile, the background flow be- 
comes super-sonic for L,- > 4/3H. This implies that for increas- 
ingly larger boxes the numerical time-step, defined by the Courant 
condition, becomes dominated by the unperturbed shear profile. 
To circumvent this undesirable constraint, it becomes mandatory 
to split-ofT the transport term due to the background profile. 

The shear transport is usually implemented in terms of an 
interpolation step. This wa s first in t roduc ed in cylindrical geom- 
etry and termed FARGO bv lMassetl l l2000h. Later, the meth od was 
adopted to the local shearing box model bv lGammid ( 120011) . 

A rather intricate extension for the induction equation that re- 
quires mapping of the magnetic field components has been pro- 
posed by iJohnson et al.l ( 120081) . We here follo w the (much sim- 
pler) c onstrained transport approach proposed bv lStone & Gardinen 
( I2OI0I) , which by construction preserves the solenoidal constraint. 

For our implementation of the orbital advection scheme, we 
operator-split the advection step from the Runge-Kutta time inte- 
gration of the remaining terms. For the interpolation of the fluid 
variables, we make use of the high-order Fourie r scheme (SAFI) 
described in appendix B of ljohansen et alj ( 120091) . 

We similarly apply SAFI to obtain the non-integer part of 
the line-integrals which contribute the circulation of the electric 
fields entering the indu ction equation (cf. Eqs. (61) and (62) in 
IStone & Gardinen 120100 . The treatment of the magnetic source 
term in the total energy equation can be successfully avoided if the 
magnetic energy is removed from the total energy during the inter- 
polation. Th e implementation has bee n tested with the advection of 
a field l oop 1 Gardiner & Stone 20081) . and the exact wave solution 
given in lBalbus & HawlevI ( l20m 

Using SAFI rather than slope-limited linear interpolation, effi- 
ciently reduces the dissipation due to the transport step, an d more- 
over its dependence on position (see ljohansen et alj|2009l) . In fact 
the scheme adds so little dissipation that the TVD requirement 
might be violated. To formally make the interpolation total vari- 
ation diminishing, we therefore discard the Fourier mode corre- 
sponding to the Nyquist frequency. 



2.1.2 Particle dynamics 

In this paper, we restrict ourselves to the study of how disc turbu- 
lence affects embedded particle populations. Neglecting their back 
reaction on the flow, particles are hence treated as passive test bod- 
ies which do not interact mutually either through physical collisions 
or gravity. Under these assumptions, we ignore the possibility of in- 
creasing the velocity dispersion of particles via mutual gravitational 
scattering. While this effect might become important for ~ 10^ km- 
sized objects, it simply adds to the external stirring. Physical col- 
lisions between planetesimals, however , can provide a source of 
damping. This effect was considered bv llda et al.l ( l2008h . and was 
found to be important in determining the equilibrium velocity dis- 
persion only for bodies with size < 1 km, and so we do not consider 
this effect in this paper. Moreover, because the particles cannot ex- 
ert drag forces on the gas, our a pproach excludes collectiv e effects 
such as the streaming instability dYoudin & Goodmanl2005i) , which 



is a focus of current numerical studie s 1 Youdin & JphansenluOOTI ; 
iBalsara et al . 2009 ; B ai & Stondl2010l : lMiniatill201(il) ^ 

We include different species of particles to quantify various 
aspects of the flow. Firstly, massless tracer particles (which instan- 
taneously follow the gas velocity) measure the Lagrangian diffu- 
sion of the flow. This is relevant for small dust grains which are 
tightly coupled to the gas. Secondly, we include particles repre- 
senting planetesimals. These particles interact with the flow via the 
gravitational potential produced by the gas density, and through the 
aerodynamic drag force. The relative importance of these effects 
is expected to change for planetesimals in the size range 1 m to 
10 km, which are the subject of this study. Finally, for the purpose 
of comparison, and as a proxy for larger objects (e.g. small pro- 
toplanets), we include swarms of particles which experience gas 
gravity but are not subject to an aerodynamic drag force. With the 
exception of the tracers, all particles are subject to the local dynam- 
ics, i.e., they experience the Coriolis force in the rotating frame and 
the tidal force stemming from the local expansion of the Keplerian 
rotation profile. As a consequence, planetesimals generally perform 
epicyclic oscillations of fluctuating amplitude, around a stochasti- 
cally migrating guiding centre. 

2.2 Numerical method - global model 

In the global simulations we solve essentially the same set of equa- 
tions for ideal MHD as described in Sect. 12.11 for the shearing 
box runs, except that we adopt cylindrical coordinates (r, (p, z) 
(see lNelsonI ( 12005 ) for a full description). The simulations are per- 
formed in a rotating frame with angular frequency equal to the Ke- 
plerian frequency at the midpoint of the radial computational do- 
main. We use a locally isothermal equation of state 



P(r) = cirY g. 



(2) 



where Cs(r) denotes the sound speed which is specified as a fixed 
function of r. T he models inve stigated may be described as cylin- 
drical discs (e.g. lHawlevl200ll) . in which the gravitational potential 
is taken to depend on r alone. Thus the cylindrical disc models do 
not include a full treatment of the disc vertical structure. Models 
of this type are employed due to the high computational overhead 
that would be required to resolve fully the disc vertical structure of 
a stratified model. 

The global simulations presented in this paper were per- 
formed using an older version of nirvana, which uses an algo- 
rithm very similar to the zeus code to solve the equatio ns of 
ideal MHD dZiegler & York3ll997l : IStone & Normadll992h . This 
scheme uses operator splitting, dividing the governing equations 
into source terms and transport terms. Advection is performed us - 
ing the second-order monotonic transport scheme jvan LeenI 19770 . 
and the magnetic field is ev olved using the Method of Characteris- 
tics Constrained Transport ( IHawlev & Stondl 19950 . 

2.2.1 Gas disc model 

The main aim of this paper is to examine the orbital evolution 
of planetesimals and smaller bodies (boulders) in turbulent discs, 
where the disc turbulence has achieved a well-defined steady state. 
Another aim is to compare local shearing box runs with global sim- 
ulations. To achieve these aims, most of the global disc models 
were chosen to have a relatively narrow radial extent, and azimuthal 
domains running between < < n/2. The turbulent stresses gen- 
erated within a global disc lead to significant changes in the ra- 
dial density distribution, such that the properties of the disc deviate 
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substanti ally from a statistical steady state over runs times of 100s 
of orbits jPapaloizou & NelsonI 120031) . In order to overcome this, 
we have introduced an additional equation to be solved alongside 
Eqs. O in the global simulations, whose purpose is to maintain a 
roughly constant surface density profile during the runs: 



dQJt) _ g(t) - go 
dt T„ 



(3) 



Here g{t) is the density at each spatial location in the disc at time t, 
go is the initial density defined at each location in the disc, and T,n is 
the characteristic time on which the local perturbed density evolves 
back toward its original value. Clearly if Tp, is shorter then any 
relevant local dynamical time in the disc, then perturbations will 
be damped very quickly and the disc will not be able to achieve a 
turbulent state. Similarly, if T,n is much longer than the longest (vis- 
cous) evolution time in the disc, the global density profile will be 
able to evolve to be very different from the initial one. By choosing 
Tn, to be intermediate between these two extremes, we find that our 
global disc models are able to develop a well defined turbulent state, 
in which the turbulence is able to maintain an approximate statisti- 
cal steady state over long time scales (> 10^ orbits). We have found 
that r„, = 50 orbits measured at the radial midpoint of the disc 
model provides a model with the desired properties. 

A possible alternative to solving Eq. ([3} would be to feed in 
mass at the outer radial boundary at the requisite rate. Such an ap- 
proach works well as a means of generating a steady dis c when 
the di sc is laminar and employs the a-model for viscosity JMassell 
l2002 n. In a turbulent disc, where there is substantial temporal and 
spatial variation in the viscous stress, however, such an approach 
may not be effective at maintaining a well-defined surface density 
profile. 



2.2.2 Initial and boundary conditions 

Each global model has a value for the inner and outer radii of the 
computational domain, ri„ and r^at, a value for the height of the up- 
per and lower "surfaces'" of the disc, Z^in and Z^ax. and minimum 
and maximum values of the azimuthal angle, 0niin and ^niax- These 
values are tabulated in table[T] along with some of the physical pa- 
rameters described below. The resolution used in each simulation 
covering nil in azimuth was (A',-, A';^, N-) = (160, 320, 40). For the 
model covering In in azimuth A'^ = 1280. This corresponds to 10 
cells per mean scale height in the radial and azimuthal directions, 
and 20 cells per scale height in the vertical domain. When consider- 
ing the behaviour of the correlation time for the fluctuating torques 
in Sect. 13.3.11 we also ran models with double and quadruple the 
resolution in the radial and azimuthal directions. 

The unit of length in our simulations is taken to be 2 au, such 
that a radial distance of r = 2.5 in computational units corresponds 
to 5 au. The unit of mass is assumed to be the Solar mass. When 
we discuss the temporal evolution of our models later in the paper, 
we adopt a time unit which is equal to the orbital period at r = 2.5 
(equivalent to 5 au), this being the midpoint of our radial domain. 

The disc models we adopt are similar locally to the minimum 
mass solar nebula model (MMSN, Hayashi 1981) at a distance of 
5 au from the star. All models have a constant aspect ratio H/r, and 
all but one model has H/r = 0.05. Model G5 has H/r = 0.075 
for the purpose of comparing directly with the local shearing box 
models (where this value of H/r was used). Initially the density is 
constant, such that the surface density S = 150 gem"'. 

The initial magnetic field in the global runs is purely toroidal, 
with the local field strength being determined from the local plasma 



P parameter (p = P^as/Pmiis)- The value of yS used in each run is 
tabulated in table [T] The field is introduced at all locations in the 
disc, except near the radial boundaries where the field is set to zero 
for r - Hn < 0. 1 and rom - r < 0.1. The initial disc velocity is 
determined according to 



GM. 



m 



(4) 



where My, is the mass of the central star, with the other velocity 
components being set to zero. Prior to a run being initiated, how- 
ever, all velocity components are seeded with random noise with 
an amplitude equal to 5 % of the local sound speed. 

Each disc model described in table[T]was evolved until the tur- 
bulence reached a quasi-steady state before the planetesimals were 
inserted. Most runs included planetesimals of size 10 m, 100 m, 
1 km and 10 km, with each size being represented by 25 particles. 
The planetesimals were distributed randomly in a narrow annulus 
centred on the computational radius r = 2.5 (equivalent to 5 au in 
physical units) of width Ar = 0.2, with their initial velocities cal- 
culated such that they are on circular orbits under the influence of 
the instantaneous gravitational potential of the central star and tur- 
bulent disc. 

We adopt periodic boundary conditions at the vertical and az- 
imuthal boundaries, and reflecting boundary conditions at the radial 
boundaries. Furthermore, we use wave damping boundary condi- 
tions in th e vicinity of the radial bo undaries using the scheme de- 
scribed in lde Val-Borro et al.l ( 120060 . This scheme relaxes the ve- 
locity and density near the boundaries toward their initial values on 
a time scale equal to 10 % of the local orbital period. 



2.3 Gravitational forces 

Both the local and the global models calculate the gravitational po- 
tential of the gas at the position of the particles via direct sum- 
mation. This is computationally favourable as long as one is only 
interested in a relatively small number of particles. Since we do 
not consider collective effects, an ensemble of a few tens of mem- 
bers is usually sufficient to reasonably determine the time-averaged 
distribution of their positions and velocities. 

For the integration of the gravitational force, the mass con- 
tained within a grid cell is treated as a point source located at 
the cell centre. To avoid artifacts due to close encounters, we 
apply a common smoothing length formalism with a parameter 
b = (6^^ + dl)'^- equal to the dia gonal across the cell. 

As can be seen in Fig. 2 of iHeinemann & Papaloizoul ( l2009bl) . 
the turbulence driven by long- wavelength MRI modes leads to the 
formation of strong density waves, which develop very little struc- 
ture along the vertical direction. Accordingly, we neglect the ver- 
tical component of the forces and adopt a cylindrical description 
where the gravity now only depends on the vertically integrated 
column density. This approach greatly reduces the computational 
demand and is consistent with neglecting the vertical density strati- 
fication. As we will see in Sect. |3.3l this 2D treatment enhances the 
gravity forces by a factor of roughly two. 

When considering local shearing box simulations in particu- 
lar, a question that needs to be answered is how large a shearing 
box is required before the gravitational forces from a turbulent disc 
are converged. In Fig. [T] we plot the relative error of the gravity 
force when integrating over spheres with increasing radius. We see 
that, for values accurate at the percent level, one requires a sphere 
of influence with a radius of about ten pressure scale heights. As ex- 
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Figure 1. Convergence study regarding the computation of tlie gravitational 
force via direct summation; the convergence order is approximately 3/2. 



pected, forces in the azimuthal direction are affected more strongly 
by long-range contributions. Because the net effect is determined 
by density fluctuations, the convergence is weaker than the r^^ de- 
pendence for Newtonian gravity. We consequently find a smaller 
convergence order of about 3/2. In our local simulations, we adopt 
a sphere (cylinder) with rd = 16//, within which we compute the 
gravitational acceleration on a planetesimal, and suitably extend 
the domain by mirroring ghost domains according to the sheared 
periodicity. 

In the global simulations the gravitational field is computed 
by summing over all grid cells. Most simulations were run using 
an azimuthal domain ofn/2, although the planetesimal orbits cover 
the full 2/r domain. When calculating the gravitational field expe- 
rienced by the planetesimals, additional copies of the disc, each 
shifted by n;r/2, where n 6 j 1, 2, 3), are used to mimic a disc which 
covers the full 2;r in azimuth. One run with a full 2;r azimuthal 
domain was run to check that the above procedure gives accurate 
results. 



2.4 Gas drag 

For the gas drag, we use the usual formulae for the Stokes and Ep- 
stein regimes ( Weidenschilling 1977; R afikov, 200 4.). In the Epstein 
regime, the drag force is given by 



with the stopping time 



Ts 



(Vg - Vp) Tj 



(5) 



(6) 



where R^, is the physical radius of the planetesimal, g is the gas 
density at the position of the planetesimal, g^, is the internal density 
of the planetesimal, Vp is the planetesimal velocity, Vg is the gas 
velocity. In the Stokes regime it may be written 



■^ drag ■ 



1 

-CunRZQ |Vp - Vgl (Vp - Vg) 



where Co is the drag coefficient, which takes the values 

24. •/?;' He < 1 

24. Kf^ !<??,< 800 

0.44 n, > 800 



(7) 



(8) 



where He is the Reynolds number of the flow around the planetesi- 
mals, defined by 



where Vpg = |Vp- Vg|, and the molecular viscosity of the flow around 
the planetesimal is given by 



= AcJ3 



(10) 



in most of the global simulations we have performed. It should be 
noted that in the local simulations, however, the molecular viscosity 
was defined by 



^cJ2, 



(11) 



so we have run one global model (G5) adopting this value. The 
molecular mean free path A = (^cth,)"', where n = Ql{p.mn) is the 
number density of particles, /j is the mean molecular weight, and 
niH is the mass of the hydrogen atom. All but one global simulation 
adopted the assumption that the disc gas is composed entirely of 
molecular hydrogen, giving // = 2. The local shearing box models 
adopted a value of /j = 2.4, so we have run one global model (G5) 
using this value in order to provide a direct comparison between 
local and global models. We adopt a value of cthj = 10~'^ cm^ fo r 
the collision cross section of molecular hydrogen jRafiko\ll2004h . 
Trilinear interpolation is used to obtain the gas density and velocity 
at the position of each planetesimal. 

The time integration of the particle motion in the shearing box 
simulations was performed as follows. Because the gas drag rela- 
tions are of the form v" in the relative velocity v = |i'g - Vp|, the 
update can be performed analytically. In the case a = 1, the rele- 
vant time scale ai"' is just the classical stopping time, and the decay 
is exponential. The update from time /„ to time f„+i with timestep 
6t can then be written as 

v„+i = v„e'"^' . (12) 

For a i^ 1 , we obtain power-law solutions for the damping of the 
relative velocity, which can be expressed as 



'{\+(a-\:)(h6t)\ 



(13) 



f^e - 2/?pVpg/yp 



(9) 



where w"' is now a generalised "stopping time" defined hy cb = 
/di 1' 'i with /dr the specific drag force acting on the particle. 

In the global simulations, t he particles were evolved using a 
fifth-order Runge-Kutta scheme jPress et al.lll996l) . 



3 RESULTS 

We organise the discussion of our results by first describing the 
evolution of the disc models. We then describe the evolution of 
the velocity dispersion (or equivalently eccentricity) of embedded 
planetesimals, followed by a discussion of their migration through 
changes in the semimajor axes. 



3.1 Hydromagnetic turbulence - local model 

Restricting one's consideration to a local approach, one might 
naively think that a relativ ely small b ox should suffice to capture 
the relevant dynamics (see lRegev & U murhan 2008, for a general 
discussion on limitations of the local approximation). As can be 
seen in Fig |2l this notion is supported by looking at the turbu- 
lent stresses created by the saturated MHD turbulence. Taken as 
the only criterion, one arrives at the conclusion that box sizes of 
<: 4// are sufficient to study the local stirring of particles - but this 
is misleading as the gravitational torques acting on the planetesi- 
mals are strongly affected by spiral density waves, which arise as a 
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Figure 2. Box size-dependence of key indicators ciiaracterising tlie liydro- 
magnetic turbulence: Wliile both the Reynolds and Maxwell contributions 
to the turbulent stress are well converged for azimuthal box sizes above 
Ly S AH. the effected density waves are clearly suppressed in too small 
boxes. Converged values for the relative density fluctuation can only be 
obtained for Ly S 8//. This is directly reflected in the gravitational torques 
acting on the particles (cf. Fig. [8). 




-2 2 
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Figure 3. Two-dimensional autocon'elation function of the vertically inte- 
grated gas density, computed from six uncorrelated snapshots of a simula- 
tion with box size 4Hxi6Hx2H. The thick black contour indicates a value 
of e"' while the thick grey line marks the first zero-crossing. 



secondary feature of the vigorously driven turbulence in a shearing 
background. 

In a series of two papers. iHeinemann & Papaloizoul ( l2009al lbr) 
study the mechanisms by which spiral density waves are excited 
in a differentially rotating fluid. One important conclusion of their 
work is that a minimum azimuthal extent L,, ^ 6H is required to 
properly capture the dynamics of spiral waves. This is exactly what 
we see when looking at the autocorrelation function (ACF) of the 
vertically integrated density, as plotted in Fig. [3] Here the thick 
black line indicates the contour where the correlation falls off to 
e"', and we see that density structures are predominantly trailing 
waves with an azimuthal extent of about six pressure scale heights. 
The aspect ratio of the approximate ellipse is about one sixth, indi- 
cating that convergence in the radial direction should be obtained 
much earlier (cf. right panel of Fig. [8}. 

Returning to Fig.|2l we see that this requirement is reflected in 
the relative rms density fluctuations ( Sg/g > measured in the satu- 
rated turbulent state (uppermost curve in Fig.|2}- Increasing the box 
size from 2H to 8//, results in an increase of ~ 75% in the relative 
rms fluctuation s. This is consider ably larger than the ~ 25% in- 
crease found bv lYang et al.l J2009I) . who performed a similar study. 
As the authors themselves mention, even this not particularly dra- 
matic effect seems to strongly affect particle stirring. This can be 
understood when taking into account the central finding of the sim- 
ulations undertaken in Hei nemann & Papaloizou (2009b), namely 
that spiral density waves quickly grow into the non-linear regime 
where they develop steep shock-like features. It seems natural that 
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Figure 4. Contributions of the Reynolds (lower) and Maxwell (middle line) 
stresses to the effective a parameter (upper line) versus time for model Gl. 



the resulting intermittent density structure creates highly fluctuat- 
ing torques, gravitationally enhancing the turbulent velocity dis- 
persion of embedded planetesimals. We examine this issue in more 
detail in Sect.|3.3|below. 



3.2 Hydromagnetic turbulence - global models 

Prior to inserting the planetesimals in the global disc models, we 
allow the MRI to develop into quasi-steady non linear MHD tur- 
bulence, with well defined volume averaged stresses operating. An 
example of the time history of the Maxwell, Reynolds and total 
volume averaged stress (taken from model GO/Gl) is presented in 
Fig. [4] and the total stress generated in each of the models as a 
function of time is shown in Fig.|5] It is clear that the models have 
evolved to a quasi-steady state. The turbulence generates a distri- 
bution of density and surface density fluctuations which are well 
fitted by Gaussian distributions. The computed standard deviations 
for these Gaussian fits (which we denote as (Sg/g) and ((5S/S», are 
tabulated in table [T] for each disc model, along with the time aver- 
aged stress parameter a. As can be seen from table [T] our models 
generate a values in the range 0.017 ^ a < 0. 101, with correspond- 
ing rms values for (Sg/g) in the range 0.101 < (Sg/g) < 0.266. 

In terms of physical parameters, the global model which is 
most similar to the shearing box simulation is run G5. We see 
that a - 0.05 in Fig. |2] for the shearing box simulation, whereas 
a = 0.035 for model G5 because of the differing magnetic field 
topologies and strengths. The density fluctuations for the shearing 
box model give {5glg) = 0.17, whereas {6g/g} = 0.13 for model 
G5. 

In Sect. l3.1l it was shown that the two-point correlation func- 
tion for the density field obtained in local shearing box simulations 
was highly anisotropic, with structure in the azimuthal direction be- 
ing stretched by a factor of ~ 6 relative to the radial direction. We 
plot the corresponding two-point correlation function for model Gl 
in Fig. [6] It is very similar to that shown in Fig.[3]for the shearing 
box run, indicating strong similarities in the density structures ob- 
tained in local and global simulations. 

3.3 Gravitational torques versus shearing box size 

To illustrate the strong dependence of the disc gravity on the do- 
main size in local simulations, we performed sets of runs with 
varying radial and azimuthal extents of the box. As a measure of 
the magnitude of the stirring, we record time series of gravitational 
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Table 1. Global simulation parameters and results. 



Simulation '■in/''out 



Zmin/Zmax HIr /3 (a) (SqIq) {61.11.) C^(iv)/[c, X 10-3] C<,(A«)/[10-''] 



GO/Gl 


1.5/3.5 


;r/2 


±0.125 


0.05 


G2 


1.5/3.5 


In 


±0.125 


0.05 


G3 


1.5/3.5 


nil 


±0.125 


0.05 


G4 


1.5/3.5 


nil 


±0.125 


0.05 


G5 


1.0/4.0 


nil 


±0.187 


0.075 



50 0.035 0.143 0.109 8.87 ± 0.50 7.92 ± 0.69 

50 0.040 0.159 0.124 8.20 ±0.32 12.20 ±0.56 

200 0.017 0.101 0.088 8.18 ±0.41 7.17 ±0.39 

12.5 0.105 0.266 0.150 11.30 ±0.70 10.00 ±0.10 

50 0.034 0.126 0.094 4.69 ±0.21 7.65 ± 0.67 




200 



Figure 5. Evolution for the total a value for the models Gl, G3, G4 and G5 
described in table [T] 
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Figure 8. Left: Specific gravitational torques vs. box size as determined 
by fitting a normal distribution (cf. Fig.|2). Gravitational torques computed 
via the column density ("2D force") are consistently enhanced by a factor 
of a roughly two, as expected from a simple geometric argument. Right: 
Isolation of the eft'ects due to variation in L^, and Lj,, respectively. 
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Figure 6. Contours of the two-point correlation of the surface density av- 
eraged over six snapshot from model Gl. The heavy contours represent the 
«- , and the zero level is represented by the grey line. 
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Figure 7. Time averaged distribution of the gravitational torques acting on 
a set of particles. Left: small box, the distribution shows excess of both 
large and small values (as typical for intermittency). Right: large box, the 
histogram is well represented by a normal distribution centred around zero 
and with a standard deviation of 1.49 X 10 cm s" . 



torques at fixed positions and compute the width of the resulting 
distribution function. This is exemplified in Fig. [T] where we see 
that, only for large-enough boxes, the torques are consistent with a 
normal distribution. Gaussian fluctuati ons, in turn , warrant stochas- 
tic modelling as considered by Youdin & Lithwick ( 2007) , for ex- 
ample, for the case of interactions via gas drag. 

The intermittent distribution seen in the left panel of Fig.|7]is 
likely related to recurring channel modes and the spiky nature of 
time series for the turbulent stresses, which occurs when going to 
too narrow boxes in the radial direction. This phenomenon, which 
is related to the truncation of the dominant parasitic modes, was 
first observed bv lBodo et al.l ( 1200 a) , and we can confirm this result 
with our simulations. 

As can be seen in the left panel of Fig. [8] the torques show 
a pronounced dependence on the box size, spanning almost one 
order of magnitude. Even for L S 8//, a weak trend towards higher 
torques is visible, albeit remaining within the error bounds. 

As already discussed in Sect. l2.3l using the column density for 
computing forces consistently enhances the torques by a constant 
geometric factor. For a fixed vertical extent L-, the results can be 
corrected accordingly. 

As expected from the shape of the autocorrelation function 
in Fig. [3] and the arguments given in Heinemann & Papaloizojj 
( l2009ar) . the observed dependence on the box size is primarily a 
dependence on L,,. This is illustrated in the right panel of Fig. [8] 
where we plot the respective dependence on Lj and Ly separately. 



3.3.1 Torque correlation time 

In addition to examining the dependence of gravitational forces on 
the box size, it is also important to consider the influence of the box 
size on the correlation time associated with the stochastic gravita- 
tional forces experienced by embedded bodies. A key issue here is 
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whether the periodicity of the shearing box, which allows for the 
possibility of waves in the flow to propagate radially past an em- 
bedded planetesimal on multiple occasions prior to damping, com- 
bined with advection due to the background shear flow, can modify 
the recurrence time of the temporally varying gravitational field. 

In the previous subsection, we discussed the value of the stan- 
dard deviation of the stochastic torques obtained for runs with dif- 
ferent box sizes. Using the time series for the stochastic torques 
obtained from these simulations, we can calculate the autocorrela- 
tion function (ACF) for the torques as a function of box parameters. 
A selection of the ACFs we obtained are plotted in Fig.[9l 

In the left panel of Fig.[9l we show the torque ACF for a box 
with size 2H x 8// in the horizontal direction. The plot can be di- 
rectly compared to Fig. 15 in Yang et al. (2009), and the results are 
very similar. The position of the first and second zero-crossing as 
well as the relati ve amplitude of the negative part of the ACF are in 
good agreement. lYang et alj have speculated that the observed un- 
dershoot may be responsible for reducing the diffusion coefficient. 
While this might in fact be the case, we suspect that the observed 
feature is probably an artifact due to the periodic boundary condi- 
tions, limited box size and aspect ratio. 

The sinusoidal modulation of the signal in the ACF is consis- 
tent with periodicity being introduced into the temporal evolution 
of the torques caused by density waves traversing the box in the ra- 
dial direction multiple times before they are dissipated. Looking at 
animations of the column density, one can identify local density en- 
hancements whenever two waves cross each other. These features, 
which create strong torques locally, are more pronounced in smaller 
boxes, and also at higher resolution (where the dissipation time is 
somewhat enhanced). Studying a set of simulations, we found the 
dependence on the box size and aspect ratio to be dominant, while 
the trend with resolution is rather weak (~ 30% when increasing 
the resolution by a factor of 3) - we hence focus on this issue. To 
pursue a more quantitative analysis, we fit the computed ACFs with 
the following model 



S r(T) = 1(1 - a) + a cos(2;r a) r)] e 



-t/tc 



(14) 



with three free parameters, namely a, indicating the relative 
strength of the proposed sinusoidal feature, a> giving its period, and 
finally the correlation time t^, assumed to be common between the 
two components. 

For the sake of simplicity, we assume that there exists only 
one wave-like mode and both components of the mix decay with 
the same characteristic time. While more elaborate fits (e.g. with 
separate decay times) produce slightly better agreement, we found 
no systematic trend in this dependence. We therefore refrain from 
the associated further complication. 

As one can see in the left and centre panels of Fig.|9l the model 
produces reasonable fits for the given data. For the small and large 
boxes, we infer relative amplitudes of the cosine-like feature of ~ 
90% and ~ 25%, respectively. This illustrates the trend with box 
size at a fixed aspect ratio. Contrary to our own expectation, the 
periodicity is not simply related to the radial extent of the box. 
This has been found studying a set of simulations, keeping L,. = 8// 
fixed and progressively increasing L,. Even at L, = 16/f, the ACF 
remains unchanged and is very similar to the one in the left panel of 
Fig.|9] It appears instead that changing both the box size and aspect 
ratio leads to changes in the torque correlation time, as illustrated in 
the centre panel of Fig.|9]for a run with box dimensions 4Hx l6Hx 
2H. The implication appears to be that the periodicity introduced 
into the run of torques versus time, and hence into the ACF, results 



from a combination of wave propagation in radius and advection in 
the azimuthal direction due to the background shear. 

Having identified a periodic feature, we can now make an 
unbiased estimation of the temporal correlation of the fluctuating 
torques. Although unexpected from a naive by-eye inspection, the 
correlation times in the two cases only differ by about 50 per cent. 
This shows that estimating the correlation time according to loca- 
tion of the first zero crossing when the autocorrelation function has 
a significant sinusoidal component can be misleading. We suggest 
that as an alternative the correlation time be estimated using a fit- 
ting formula such as given in Eq. ( 114b . We see from the central 
panel of Fig.[9]that r^ = 0.32 for our large shearing box model. 

Finally, in the right panel of Fig. [9l we compare the torque 
ACFs of our global simulation Gl with local simulations at compa- 
rable resolution. We see that excellent agreement is obtained when 
using large-enough boxes with an elongated aspect ratio, indicat- 
ing that the correlation time for the global simulation Gl t^ - 0.32. 
Interestingly, the correlation time measured from the zero crossing 
point for model Gl is t^. = 0.68, about a factor of two larger than 
we infer from the fitting procedure described above. 

In order to check the sensitivity of model Gl to numerical res- 
olution, we have re-run it at both double and quadruple resolution 
in the radial and azimuthal directions, giving (Nr x A'^) equivalent 
to (320 X 2560) and (640 x 5120), respectively for a disc which 
covers the full 2n in azimuth. These runs have 20 and 40 cells per 
mean scale height in the radial and azimuthal directions. The ACFs 
measured in each of these runs are very similar to that shown in the 
right panel of Fig.|9l 

3.4 Gravitational stirring versus distance in global models 

As described in Sect. 13.21 the disc models listed in table [T] were 
evolved until they had reached a statistical steady state, prior to in- 
serting the planetesimals. Before discussing the evolution of plan- 
etesimals which experience gas drag, we examine the length scales 
over which gravitational stirring of planetesimals by the turbulent 
disc occurs. 

We performed a series of simulations using the disc model 
GO, where we defined a gravitational sphere of influence of varying 
size, Rcui, around the planetesimals (which do not experience the 
gas drag force in this particular simulation) . We then examined 
how the radius of this sphere of influence changed the evolution 
of the planetesimal velocity dispersion. The sphere of influence in 
each simulation is defined by a radius, /?cut, measured in units of 
the local value of H. Gas within this sphere of influence exerts a 
gravitational force on the planetesimal. At the edge of this sphere, 
the contribution to the force is tapered to zero over a distance equal 
to H using the hyperbolic tangent function. For example, a sphere 
of influence equal to iJ^ut = 2// allows a full contribution to the disc 
force within IH, and beyond this the force contribution is tapered 
to zero. 

Fig. [To] shows the evolution of the rms of the radial velocity 
dispersion, cr(vr), as a function of ifcut, where (r{v,) is measured in 
units of the sound speed. For this model, c^ = 666 m s"' at r = 2.5 
(5 au). It is clear that contributions to the gravitational force occur 
even beyond a cut-off radius R^ut = 8//, in agreement with the cal- 
culations presented in Sect. |3.3| for shearing boxes of different size. 
The basic reason for this has already been explained: spiral waves 
generated by the turbulence provide coherent structures which are 
stretched in the azimuthal direction by the shear, and contribute sig- 
nificantly to the stochastic gravitational forcing experienced by the 
planetesimals. Clearly, computational domains are required which 
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Figure 9. Box size-dependence of the torque autocorrelation. Left 6- Centre: Best fits according to Eq. J 141 (solid black line) for a small and large box, 
respectively. The torque ACFs (dark grey lines) are measured at a fixed position, and error intervals (shaded areas) are estimated from considering sub- 
intervals in time. The first and second zero-crossing ai'e indicated by triangles. Right: Comparison between local models and the global ran Gl applying 
comparable spatial resolution. Excellent agreement is obtained if a large enough box-size is chosen. 
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Figure 10. Evolution of the radial velocity dispersion, (t(iv), in units of the 
local sound speed, as a function of the size of the gravitational sphere of 
influence described in the text. 
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Figure 11. Exemplified temporal evolution of the radial displacement Ax 
for a swarm of eight gravitationally excited particles. The initial positions 
have been off-set for clarity as indicated by the dashed lines. 



are large enough to capture the gravitational stirring which occurs 
on scales up to 8 - 10 scale heights, at least in the azimuthal direc- 
tion. 

Examining the curve in Fig. [TO] which corresponds to no cut- 
off in the gravitational force, we see that the evolution of cr(vv) can 
be reasonably well fitted as a random walk. The smooth solid line, 
and the dashed lines, correspond to the function Co- yjt - to, with 
Co- = 8 ± 0.8 X 10"\ where the time is measured in orbits at r = 2.5 
(5au). 



3.5 Evolution of the planetesimal velocity dispersion - local 
model 

After to = 20 orbits, when the turbulence driven by the magneto- 
rotational instability has reached a quasi steady-state, we disperse 
several swarms of test particles into the flow. For easy reference, we 
label these sets as "G" for the particles that experience the gas grav- 
ity only, "D" for particles subject to the gas drag-force, "G+D" for 



the combined effect, and "T" for massless tracer particles. While 
the sets G and T consist of only one species, the sets D and G+D 
are composed of ten species each. Particle radii are R = 1 m, 2 m, 
5 m, 10 m, 20 m, 50 m, 0.1km, 0.2 km, 0.5 km, and 1km, respec- 
tively. Eight particles are used to represent each size. Due to the 
low number of particles, we expect sampling errors on the order of 
20 - 35%. In reality, these numbers have to be seen as upper limits. 
Because our fits are based on time histories, the effective statisti- 
cal basis is probably somewhat larger than for any given instant 
in time. To roughly quantify the uncertainty due to Poisson fluctua- 
tions, we performed a lower resolved fiducial run with 100 particles 
for the gravity-only set. Taking twelve subsets of 8 particles each, 
we arrive at a standard deviation of 23% amongst the different re- 
alisations. For six sets of 16 particles, this number is only slightly 
reduced to 18%, justifying the initial choice. 

Depending on the particle size and the prevailing form of the 
coupling, the time evolution of the various species is quite diverse. 
Massless tracers and small particles with R^, ^ 10 m, whose dynam- 
ics are largely controlled by gas drag, essentially follow the turbu- 
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Figure 12. Illustration of the algorithm used to measure the eccentricity 
{lower panel) of a particle moving on epicyclic orbits via tracking the aphe- 
lion and perihelion {upper panel). 
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Figure 13. Evolution of the radial velocity dispersion for set 'G+D'. For 
sizes above S i 1 km, forces due to gas drag become gradually negligible 
and the evolution is described by a V? behaviour 



lent flow and describe a random walk. The larger G+D particles, 
as wefl as the G set are coupled more weakly. For them the orbital 
dynamics due to the Coriolis and tidal forces becomes increasingly 
relevant. Viewed from a local perspective, the motion of these par- 
ticles can be described as epicyclic oscillations with a modulated 
amplitude and a stochastically migrating guiding centre, as seen in 
Fig. \TT\ Thus, the relevant properties of the particles' motion can 
be described using two characteristic quantities: (i) orbital eccen- 
tricity, or the amplitude of the epicyclic motion (or equivalently the 
velocity dispersion relative to a circular Keplerian orbit); (ii) semi- 
major axis ~ the position of the guiding centre, which evolves as 
a random walk as the particles migrate from their initial locations. 
In this section we examine the eccentricity/velocity dispersion, and 
consider the migration in Sects. 1377] and |3. 8.2] 

To separate the stochastic motion of the guiding centre from 
the epicyclic oscillation, we box car-average the velocity apply- 
ing a filter-scale equivalent of the orbital frequency. We compute 
the eccentricity of the particle's orbit by tracing the position of the 
aphelion and perihelion, respectively. From these, we compute 



0.25 1 



Rq + Xp 



and 



k-l 
k+l 



(15) 



with Rq = 5 au the location of our local box. The tracking of 
the extrema and the resulting eccentricity function are illustrated 
in Fig. [T2I where we see that e(t) itself follows a random walk. 
This means that gravitational stirring can both excite and damp the 
epicyclic motions of individual particles. As we will see shortly, 
however, the influence of the stochastic forcing on the ensemble of 
planetesimals will increase the rms eccentricity amongst the mem- 
bers, and therefore heat-up the ensemble as a whole. 



3.5.1 Saturation amplitudes 

Needless to say that the velocity fluctuations will not grow ad in- 
finitum, but will reach a saturated state once the aerodynamic damp- 
ing reaches the level of the stochastic forcing. This is illustrated in 
Figs. [T3] and [141 where we plot the temporal evolution of the radial 
velocity dispersions for large and small bodies, respectively. 

Looking at the lowermost curves in Fig. [T3] we see that the 
random velocity of particles with 50 m S Rp S 200 m saturates on 
time scales of several hundred orbits. In this regime, the satura- 
tion amplitude increases with the particle size, reflecting the weaker 
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Figure 14. Same as Fig. 1131 but for smaller particles, subject only to gas 
drag. Note the opposite trend with particle radius (also cf. Fig. llSl belovy). 



damping for larger bodies. Because the gravitational stirring is in- 
dependent of the particle mass, smaller bodies reach their saturated 
state earlier compared to heavier ones. For larger bodies (500 m, 
1 km), we do not achieve equilibrium, but after 500 orbits we note 
that it(v,) ^ O.OSCs for Ikm-sized bodies (where c^ = Ikms"' in 
this model). 

Unlike for their larger counterparts, smaller bodies are pre- 
dominantly affected by gas drag and we observe the opposite trend 
with respect to the saturation amplitudes (see Fig. I14t . At first, it 
seems surprising that cr(v^) can be excited beyond the turbulent ve- 
locity Vrni.s = 0.13 Cs of the gas. This is, however, only the case for 
the radial velocity and can be understood from the orbital dynamics 
which le ads to ( v^ ) x 2{\\ ) as ch aracteristic of epicycles (also cf. 
Fig. 5 in lYoudin & Lithwickll2007l) . 

The saturation amplitudes in the velocity dispersion which 
arise as a function of particle size are compiled in Fig. [15] With 
the exception of the Rj, = Im species, which is tightly coupled to 
the gas, cr,(Vv) ~ 2tTs(v,.) as expected from orbital dynamics. The 
results for set D can be readily compared to the upper pane|j of 



' As we will discuss shortly in section |3.8.2| we find a characteristic eddy 
time Te » il"' in our simulations such that the standard case applies. 
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particle size rises in Fig. \\5\ but the largest bodies we consider 
have not had time to reach the equilibrium values of their veloc- 
ity dispersion. The explanation for the existence of the minimum 
velocity dispersion observed in Fig.[T5]is straightforward. Smaller 
bodies are tightly coupled to the gas, and so attain a velocity dis- 
persion close to that of the gas itself. Larger particles experience 
a significantly smaller gas drag, which then contributes weakly to 
counterbalancing the stirring effect of the stochastic gravitational 
force provided by the disc, leading to a large velocity dispersion. 
Intermediate sized particles in the 50 - 100 m range are sufficiently 
decoupled from the gas that they experience an orbit-averaged drag 
force that causes significant damping of the eccentricity growth 
driven by the disc gravity. 



1 '0 ' 00 

ploneteslmal radius R [m] 

Figure 15. Saturated velocity dispersions as function of planetesimal radius 
i?p, and stopping time Tj, respectively. Particles which feel both the gas 
gravity and drag force (symbols) deviate from the ones that only feel the 
gas drag (solid tines) for sizes larger than i? ai 50 m. 
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Figure 16. Evolution of the radial velocity dispersion, cr(Vr), in units of the 
sound speed, for planetesimals of different size from run G5. 



Fig. 5 in I Youdin & Lithwickl ( 120071) ■ which shows saturation am- 
plitudes obtained from a simple turbulence model. Considering the 
vast differences in the two approaches, the results agree surpris- 
ingly well. 

Considering the G-l-D set, the particle rms velocity takes its 
maximum for t^H and falls-off to about 10% its peak value 
at Ts ^ 10^ fl~'. This is a some what shallower decline than in 
the model of 1 Youdin & Lithwickl where this value is reached at 
10^ already. The minimum radial velocity dispersion obtained is 
- O.lvinis (where Vi„,s is the turbulent velocity dispersion of the gas). 
This minimum value marks the point where the G+D set deviates 
from the D set (see rhs in Fig.1151). This implies that objects of size 
»: 50 m enjoy the relative comfort of being the least affected by their 
turbulent surroundings, with a velocity dispersion that corresponds 
to 15 - 20ms"'. Interestingly this is close to the speeds required 
for catastrophic break-up of planetesimals in the 10 - 100 m range 
JBenz & Asphaug|ll999l : IStewart & Leinhardlll2009l) 

For larger sizes, the velocity dispersion is seen to rise as the 



3.6 Evolution of the planetesimal velocity dispersion - global 
model 

We now examine the evolution of the velocity dispersion of plan- 
etesimals of different size which experience the gas drag force 
within global models, focusing on the radial component of this 
quantity. The key issues that we explore are the magnitude of the 
velocity dispersion attained as a function of planetesimal size and 
disc model parameters, and the implications for the outcomes of 
collisions between planetesimals for the growth of planets. 

3. 6. 1 Model G5 - comparing local and global models 

The global disc model with physical parameters most similar to the 
shearing box simulation is model G5. As described in Sect. 12.41 
the shearing box model adopted H = 0.075, and used Eq. dill ) 
and /J = 2.4 in defining the strength of the gas drag force. Most of 
the global simulations adopted model parameters which are slightly 
different from these, a fact which was discovered after most of 
the simulations presented here had been completed. We present 
here, however, a global model with the same parameters used in 
the shearing box run for the purpose of providing a direct compari- 
son. The only difference in the underlying disc models is the choice 
of magnetic field topology and strength, with the resulting a values 
being a =^ 0.05 for the shearing box and a == 0.035 for the global 
model, which are similar enough for a meaningful comparison to 
be made. The planetesimal sizes considered in this run were 10 m, 
20 m, 50 m, 100 m, 1km and 10 km, with each size being repre- 
sented by 25 particles. 

The evolution of the radial velocity dispersion, expressed in 
units of the local sound speed, is shown in Fig.[T6]for planetesimals 
with sizes in the range 50 m - 10 km. Comparison with Fig. 1131 
which shows the same data for the shearing box simulation, indi- 
cates that there is good agreement between the two simulations. 
The 50 m and 100 m sized objects quickly attain equilibrium val- 
ues for cr(vr) in the range 0.01-0.02cs (where c^ = 1 km s"' for 
this model), similar to the values obtained in the shearing box run. 
After 500 orbits cr(vr) - O.OScj for the 1 km sized bodies, and the 
10km bodies have cr(v,) =: O.lCs, again in good agreement with the 
shearing box results. 

In principle we would expect the shearing box simulations to 
generate slightly larger velocity dispersions, due to the more vig- 
orous turbulence exhibited by that model, and the correspondingly 
larger value of (Sg/g). The fact that the gravitational stirring is ac- 
tually found to be very similar is probably because the global sim- 
ulations allow for a faster growth of the velocity dispersion due to 
larger length scales being included in the disc gravity force calcula- 
tion. More distant density perturbations are thus able to contribute 
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Figure 17. Saturated values of it(i',-), measured in metres per second, as a 
function of particle size from run G5. 



to the stochastic gravitational field experienced by the planetesi- 
mals, providing a small boost to the gravitational stirring. 

Whereas the larger planetesimals have not achieved an equi- 
librium value for it(v,-) by the end of the simulation, planetesimals 
with sizes in the range 10 m - 100 m have. We plot the saturated 
value of a(Vr) as a function of particle size in Fig.[T7l and in agree- 
ment with the results obtained for the shearing box simulation, we 
observe that there is a minimum value for cr(vv) for bodies of size 
- 50m, corresponding to a(yr) — 20 ms"'. 



3.6.2 Evolution as a function of a: models Gl - G4 

We now consider the evolution of the velocity dispersion as a func- 
tion of the turbulent strength, as measured by a, by presenting the 
results from models Gl, G2, G3 and G4. As shown in table [T] the 
value of a was modified by changing the strength of the net toroidal 
magnetic field in the initial conditions. Models Gl and G2 difi'ered 
only in the size of their azimuthal domains (;r/2 for Gl and 27r for 
G2), but the increase in domain size for G2 causes a small increase 
in a (presumably due to the presence of low-order MRI modes). 
We note that the disc parameters and values used in computing the 
gas drag strength were difi'erent in these models compared with run 
G5, and so these runs are not directly comparable with that one. 

The evolution of crCv,) for models Gl - G4 are presented in 
Fig. [TS] for planetesimals of size 10 m, 100 m, 1km and 10 km. 
Overall, the evolution of (t{\\) is found to be a weak function of 
a (we find an approximate scaling cr(\',^ oc o-"-^", see below), and 
the evolution of cr(vr) for models Gl and G2 are in good agreement. 
Although a (and {61,/1,}) is larger in model G2, we find that the ve- 
locity dispersion increases at a slightly slower rate for the larger 
planetesimals than found in model Gl. The difference, however, 
is well within the ^/N variations expected for the low numbers of 
particles used. 

In comparison with run Gl, G3 shows slower growth of cr(iv) 
for the larger planetesimals, and a smaller saturated value of (t(vv) 
for the 10 m sized bodies, as expected given the smaller value of a. 
Run G4 shows correspondingly faster growth, and larger saturated 
values, of a-(vr) due to the larger value of a. Each of the plots in 
Fig.[T8]show fits to the data for the 10 km bodies, assuming a func- 
tional form <T(Vr) = Ca-(vr) yft. The values of Ca-iv^) for each model 



are tabulated in table [T] Fitting the data for a and Ca-(vr) listed in 
table[T] using an expression of the form Ccr(Vr) = K^,^a'', leads to a 
best-fit solution with q = 0.20 and A",,, = 1.64 x 10"^. We use this 
fit to the data in Sect. 13.6.41 below, where we discuss the expected 
saturation value of (t(v, ) for 1 km and 10 km sized planetesimals as 
a function of a. 



3.6.3 Saturation values ofcr(Vr)for 1km and 10km planetesimals 

Our simulations have not run for sufficient time for (t(v, ) to reach 
its equilibrium value for the 1 km and 10 km sized bodies. Assum- 
ing that saturation is reached when the stochastic gravitational forc- 
ing is balanced by gas drag damping, we can estimate the saturation 
values by equating the forcing and damping time scales. Work- 
ing in terms of the orbital eccentricity (where e =^ Vdisp/v<., with 
I'disp = (T(vr) and vk being the Keplerian velocity), and using the 
expression Vdisp = Cv,. y/i we can write 



le^vi 



de/dt Q(v,.)2 ' 



(16) 



where T„r„„ is the eccentricity growth time. The damping time for 
the velocity dispersion can be estimated simply from the ratio of 
the momen tum associated with the velocity dispersion and the gas 
drag force llda et alj2008l) 
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(17) 



Equating expressions l |16t and l |17t , and writing the planetesimal 
mass in terms of its radius and internal density, pp, leads to the 
following expression for the equilibrium velocity dispersion 



Vdisp 



iCoQ 



(18) 



Noting that Cs = 666 ms"' at 5au in our global disc models Gl~ 
G4, and expressing Co-(Vr) in S.I. units, we obtain C,.,. = 3.15x10""* 
for model Gl. This leads to estimates of the equilibrium veloc- 
ity dispersion of I'disp = 765 ms"' for 10 km bodies (approxi- 
mately the sound speed), and Vdisp = 356 ms"' for 1km bodies 
with jpp = 3 gem"'' at 5au in a disc with E = 150 gem"- and 
Hjr = 0.05. These values are clearly very much in excess of the 
velocities required for catastrophic disruption of I km and 10 km 
sized planetesimals, being many times larger than the escape ve- 
locities, Vesc, from these bodies (vesc - 12 ms"' for a 10 km body 
with pp = 3gcm"''). Extrapolating forward in time, the time re- 
quired to reach the equilibrium value for Vdisp for the R^ = 1 km 
planetesimals is ^ 1.75 x 10^ yr, comparable to the runaway growth 
time scale at 5 au. 



3.6.4 Saturation of(T(Vr) as a function of a 



Using the expression Ccriv,) = A',,,.q-''-'' discussed in Sect. 13.6.21 in 
conjunction with Eq. l llSt , we can estimate the value of a which 
leads to a particular value of Vdisp for a particular size of planetesi- 
mal 



Vdis 



4g^R^Kla ' 



,0.4 N 1/3 



(19) 



Using the results from model Gl, and working in S.I. units, we 
obtain AT,,,. = 6. 16 x lO"**. Values of Vdisp as a function of a are 
plotted in Fig.[T9|for planetesimal sizes 1 km and 10 km, and for 
planetesimal internal densities Pp = 1 and 3 g cm"^ . It is clear that 
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Figure 18. Evolution of the radial velocity dispersion, it(v,-), in units of the local sound speed, for planetesimals of diiferent size from runs Gl (top left panel), 
G2 (top right panel), G3 (lower left panel), and G4 (lower right panel). Each panel also displays a fit to the random walk behaviour, as described in the text. 



even with a small value of a = 10"*, equilibrium velocity disper- 
sions are in the range 60 < Vdisp < 110 ms"', somewhat larger than 
the escape velocities and catastrophic disruption velocities of these 
bodies. Evidently a protoplanetary disc needs to be very quiescent 
near its midplane in order to allow for runaway growth to proceed, 
and to prevent the catastrophic breakup of colliding planetesimals 
in the 1 - 10 km size range. 

There are obviously a number of caveats contained within the 
simple arguments presented above. In a disc in which the density is 
vertically stratified (with or without a dead zone), the dependence 
of ^disp on a may be steeper if the stochastic forcing of planetes- 
imals depends on the magnitude of local density rather than sur- 
face density perturbations. It is also possible that mass accretion 
through the disc may be generate d by the winding up o f net radial 
fields in a disc with a dead zone ( Turner & Sandl200a) . breaking 
the link between a and (Sg/g), which the above arguments rely 
upon. Confirmation of the result obtained above will be sought in a 
forthcoming paper in which we examine planetesimal dynamics in 
vertically stratified discs with dead zones. 



3.7 Radial migration of planetesimals - global models 

We now consider the radial drift of the planetesimals, arising both 
from the effects of gas drag and from the action of the stochastic 
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Figure 19. Variation of the saturated velocity dispersion as a function of the 
turbulent a and the planetesimal density. 



gravitational torques. The radial drift of large planetesimals is po- 
tentially an important effect during planet formation, as it may lead 
to the delivery of material to regions which have suffered from de- 
pletion, and may also enhance the delivery of volatiles to the inner 
system from beyond the snowline. Radial drift due to gas drag is 
not present in the local shearing box simulations, since there is no 
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radial pressure gradient to generate sub-Keplerian velocities in the 
gas. We therefore begin our discussion by examining the results of 
the global simulations. 

In a sub-Keplerian disc, the gas drag force given by Eq. ([7} 
will lead to a loss of angular momentum and radial drift of planetes- 
imals at a rate which depends on their size ( Weid enschilling 19770. 
The time over which the gas drag removes an amount of specific 
angular momentum equal to A j from an embedded planetesimal is 
given by 



Tdri 



8i^p£)p Aj 
3CdQ r^\% ' 



(20) 



where g^ is the density of the planetesimal material (here assumed 
to be 3 g/cm^^). Superposed on this inward drift are the stochas- 
tic torques experienced by the planetesimals in a turbulent disc, 
which will cause diffusion of planetesimal semimajor axes (as well 
as contributing to the eccentricity evolution). The effective diffu- 
sion coefficient associated with the diffusion of planetesimal angu- 
lar momenta can be approximated by Dj = cr^Tcorr, where ar is the 
standard deviation of the distribution of torques (assuming Gaus- 
sian statistics), an d Tcorr is the correlati on time associated with the 
stochastic torques jjohnson et alj2006h . The time scale over which 
diffusion will change the specific angular momentum of a typical 
planetesimal by an amount equal to Aj is 



Tdiff ■ 



(A;y 



(21) 



We expect that the effects of radial drift and diffusion will be 
comparable when eqs. l l20t and J21t are equal, and this is the time 
scale over which the drag-induced inward drift of planetesimals of 
size /?p should just become apparent against the isotropic diffusion 
generated by the turbulence. The random walk associated with the 
stochastic torques means that the planetesimal size, i?p, for which 
inward drift just becomes apparent depends on the magnitude of 
the change of angular momentum, Aj. Larger values of Ay imply 
longer evolution times, such that the net inward drift of larger plan- 
etesimals, which are increasingly impervious to gas drag, eventu- 
ally becomes apparent. 

Results from model Gl are shown in Fig. [20] which shows 
the relative change in semimajor axis versus time for all planetes- 
imals, and it is clear that the Rp = lOm bodies undergo gas drag- 
induced inward migration, whereas the larger bodies are dominated 
by stochastic migration. The rms relative change in semimajor axis, 
a(Aala), for the R^ = 10 km planetesimals is shown in Fig.l21l 

Small changes to the specific angular momentum and semi- 
major axis of a body are related according to 



(22) 



A; _ 1 Aa 
;■ 2a' 

Using Eqs. ( 121b and ( I22t . we can estimate how far we expect large 
planetesimals to have diffused in semimajor axis during the 1200 
orbits of run time in model Gl. The distribution of torques, aver- 
aged over all 10 km-sized planetesimals, is presented in Fig. 1221 
and the standard deviation of this torque distribution is found to be 
(Tr = 1.5x 10* cm^ s"^ (equivalent to 3.46x 10"' in code units). The 
correlation time, Tcorr, was estimated to be Tcon = 0.32 orbits from 
the exponential decay of the ACF shown in Fig. [9] in Sect. |3.3.T| 
The diffusion coefficient, Dj, expressed in code units is given by 
Dj = 9.51 X 10 ', which leads to the prediction that a typical plan- 
etesimal located at r = 2.5 (5 AU) in model Gl will diffuse such 
that Aa/a = 0.021. Comparing this with the value of cr(Aa/a) in 
Fig-Hi] we see that the simulation resulted in (T(Aa/a) ^ 0.03 after 
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Figure 21. Evolution of the rms change in semimajor axis for 10 km sized 
planetesimals, where each line is annotated with the ran label. 



1200 orbits, in reasonable agreement with the value predicted by 
the diffusion coefficient. The fact that the estimated amount of dif- 
fusion is too small by a factor of ~ yl compared to the simulation 
result indicates that the approximation to the diffusion coefficient 
given by Dj = cr^T^gi-r is too small by a factor of ^ 2. This possibly 
arises because of the ambiguity discussed in Sect. 13. 3.TI regarding 
the definition of the correlation time, but may also be affected by 
sampling errors that arise from using only 25 particles per size bin, 
leading to Poisson errors at the ~ 20 % level. 

We can now examine if the gas drag-induced radial drift of the 
i?p = 10 m bodies observed in Fig.|20] and the stochastic migration 
of all larger planetesimals, is expected. According to Eq. l l22t and 
Fig. [2T] stochastic migration causes a relative change in specific 
angular momentum Aj/j = 0.015. Putting this value into Eq.OOIl 
tells us that planetesimals of size Rp ^ 25 metres should radially 
drift due to gas drag and diffuse by a similar distance. Planetesimals 
which are smaller than this should show strong inward drift, and 
larger bodies should have their evolution dominated by diffusion. 
The relative change in semimajor axes for all planetesimals in run 
Gl plotted in Fig. [20] show good agreement with this expectation. 

Changes in the mean drag-induced radial drift rate due to 
the turbulence can only be discerned for planetesimals of size 
/?p = 10 m in Fig. [20] where the straight black line plotted in the 
top left panel shows the trajectory of a planetesimal embedded in 
a laminar disc. Up to a time of approximately 300 orbits, the mean 
drift rates of the planetesimals in the turbulent disc closely match 
the trajectory in the laminar disc. At later times, however, the tra- 
jectories diverge, with the drift rates for the planetesimals in the tur- 
bulent disc decreasing by about a factor of three. This is not due to 
the effects of stochastic forces, which earlier on in the evolution are 
seen to cause a dispersion of the trajectories around a mean which 
matches the laminar case closely. Instead, variations in the effective 
a stress parameter as a function of radius cause radial structuring 
of the disc, such that peaks and troughs in the mean surface den- 
sity arise. The slowing down of the radial drift observed in Fig. [20] 
arises when the planetesimals enter a region of the disc in which 
the magnitude of the radial pressure gradient decreases. Thus, we 
observe that significant changes to the radial drift of planetesimals 
do not occur due to stochastic forcing, but may arise because of 
radial structuring of the disc due to spatial variations in turbulent 

stresses, which persist over the run times of the simulations^ 

^If we adopt a typical disc life-time of 5Myr (Hai sch et alj 

I2OOII) . then we can estimate the amount of radial diffusion for 
large (Rp > 1 km) planetesimals during the planet forming epoch. 
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Figure 20. Evolution of the semimajor axis for all planetesimals in run Gl. The top left panel displays data for the 10m planetesimals, along with the 
theoretically expected drift within a laminar disc (solid black line). The top right panel coiTesponds to the 100 m planetesimals. The lower left panel refers to 
the 1 km planetesimals, and the lower right panel shows the 1 km bodies. 
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Figure 22. Distribution of torques experienced by the 10 km sized planetes- 
imals in run Gl, averaged over all planetesimals. The standard deviation 
obtained for the distiibution o" = 1.5 X 10^ cm^ s"'. 



Adopting a diffusion coefficient, Dj, which gives agreement with 
the amount of radial diffusion observed for 10 km-sized planetes- 
imals in Fig. [21] Ea.(l2U predicts Aj/j = 0.3, corresponding to 
a 50 % change in the semimajor axis of a typical planetesimal lo- 
cated initially at 5 AU. Considering the relative contributions to 
radial drift of solid bodies from both gas drag-induced migration, 
and diffusion due to the stochastic gravitational field of the disc, we 
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Figure 23. Variation of rjrag and rjifj as a function of planetesimal size. 
The underlying disc model is the same as in run Gl . 



have plotted the evolution times Tj^ig and Tjiff in Fig. [23] assuming 
a 30 % change in the specific angular momentum of a planetesi- 
mal, using Eas.(l20t and (I2U. We can see that for evolution times 
of 5 Myr, the drag and diffusion time scales are equal for planetes- 
imal size R^ ^ 100 m, but drag time scales are much larger than 
diffusion time scales for larger bodies such as 1 km and 10 km- 
sized planetesimals. This shows that large scale migration of large 
planetesimals over evolution times of 5 Myr will be dominated by 
diffusion and not gas drag-induced inward drift. 

Strong radial mixing of planetesimals at the level of Aa/a = 
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0.5 would have very significant implications for planetary forma- 
tion. In our own Solar System, for example, radial mixing of large 
icy planetesimals from beyond the snowline would substantially in- 
crease the volatile content of the terrestrial planets, and effectively 
smear out the observed apparent variation in asteroid properties as 
a function of heliocentric distance. It thus seems unlikely that such 
radial mixing occurred in the Solar nebula, providing strong cir- 
cumstantial evidence that the degree of turbulence present in model 
GI was not present in the Solar nebula during the formation of the 
Solar System. We explore the degree of radial mixing as a func- 
tion of turbulent strength (measured by the effective a parameter) 
below. 



3. 7. 1 Planetesimal diffusion as a function of a 

The rms value of the relative change in the semimajor axes for 
10 km sized planetesimals for various models is shown in Fig.l2II 
The discs in models G 1 and G5 gave rise to very similar values of a 
and {5S/S>, and we see that the radial diffusion rates are also very 
similar. Model G3 had a value of o- = 0.017, and a correspondingly 
smaller value of ((5S/S>, which leads to a slower rate of diffusion. 
The rate of diffusion in semimajor axis is a fairly weak function of 
a (see discussion below), such that it is only after a time oft~ 600 
orbits that the divergence in the diffusion rates can be detected for 
models GI and G3 in Fig.[2T| Model G4 has a substantially larger 
value of a = 0. 1 1 , and leads to a noticeably larger rate of diffusion 
of Aa/a in Fig. 12 1 1 

We have fit the data for cr(Aa/a) shown in Fig. [21] assuming 
a simple random walk model o-(Aa/a) = Ccr(Aa) yfi, and the co- 
efficients Co-(Aa) are given in table \T\ If we fit the data given in 
table[T]for a and Co-(Aa) using the function C„-(Aa) = K/^a'', then 
the best-fit solution gives q = 0.20 and ATao = 1.70 x 10"^ This 
is the same power law dependence found for the evolution of the 
velocity dispersion as a function of a. This fit allows us to estimate 
the amount of radial diffusion that we expect as a function of a. 



3.7.2 Radial diffusion constraints on a in the solar nebula 

As discussed above, large planetesimals experience only small 
changes in their angular momenta and semimajor axes due to gas 
drag, assuming that they remain on approximately circular orbits. 
For example, a 1 km sized planetesimal orbiting at 5 au in model 
GI, in the absence of turbulence, would migrate a distance of 
~ 4 X 10"^ au in 5Myr. For such large bodies, evolution of the 
semimajor axis is dominated by turbulence-induced diffusion, es- 
pecially if the disc is fully turbulent as in models GI - G5. Model 
GI predicts that a typical 1 0km sized planetesimal orbiting at 5 au 
will diffuse a distance Ac - 2.5 au in 5 Myr. If this 50 % change 
in semimajor axis is applied to bodies in the asteroid belt, then this 
level of diffusion is probably inconsistent with the variation in ob- 
served properties of aster o ids as a function of helioc entric distance 
JGradie & Tedescolll98A iMothe-Diniz et alj|2003l) . These varia- 
tions would have been smeared out substantially if such a large 
amount of orbital diffusion had taken place during Solar System 
formation. 

The survey of how the distribution of asteroidal taxonomic 
types varies as a functio n of heliocentric distance presented by 
iGradie & Tedescd ([1982) concluded that there is a systematic vari- 
ation. For the seven taxonomic classes identified in their sample, it 
was suggested that the distribution of each peaks at a different loca- 
tion in the asteroid belt, with a dispersion around this peak location 



of 0.5 - 1 au. The observed correlation between asteroidal type 
(assumed to relate to composition) and heliocentric distance was 
interpreted as being evidence that the asteroids were formed essen- 
tially in their observed locations, being subject subsequently to a 
relatively modest amount of ra dial mixing at th e level of ~ 0.5 au. 
The more recent survey of Mothe-Diniz et alj l|2003|) presents a 
more complicated picture in which S-type asteroids are distributed 
more uniformly throughout the asteroid belt, but are nonetheless the 
dominant class in the inner and middle belt, with the more volatile- 
rich C-type bodies being preferentially located in the outer belt be- 
yond ^ 3 au. 

Radial mixing of planetesimals during planet formation is ex- 
pected to occur through gravitational interaction with planetary em- 
bryos. Wetherill ( 1992) suggested that a population of ~ Mars-mass 
embryos embedded within the primordial asteroid belt could ex- 
plain the required mass de pletion and radial mixing. More recen t 
simulations of this effect dPetit et al.l 1200 ll : lO'Brien et al] l2007h 
show that radial diffusion of asteroids peaks at a value Aa =: 0.5 au, 
broadly consistent with the observational constraints. In a series of 
related simulations, lO'Brien et al.l yOOg) examined the formation 
of terrestrial planets, and estimated the rate of water delivery to the 
Earth by accretion of volatile-rich material from beyond ~ 2.5 au. 
It was found that simulations initiated with Jupiter- and Saturn- 
analogues on circular orbits could reasonably explain the abun- 
dance of water on the Earth (^ 5 x 10"* Me). Substantial radial 
mixing of volatile-rich planetesimals and embryos from the outer 
asteroid belt, due to disc turbulence, would very likely lead to ter- 
restrial planets which are endowed with water and other volatiles 
well in excess of what is observed. 

We suggest that orbital diffusion due to turbulence which ex- 
ceeds Aa ^ 0.5 au is probably inconsistent with observations, and 
we now discuss which value of a is consistent with Aa ^ 0.5 au in 
the asteroid belt. It should be noted, however, that using the aster- 
oid belt as a test for dynamical mixing requires that the observed 
variation of propertie s with heliocentric dis tance is primordial. The 
recent suggestion by i Levison et alj ( 120090 that cometary material 
originating from further out in the Solar System became implanted 
in the asteroid belt during the late heavy bombardment means we 
should take care not to over interpret the observations, or their the- 
oretical implications. 

Using the approximate fit Cah = ^aoQ'''"'', discussed in 
Sect. 13.7.11 combined with cr{Aa/a) = C^a V?, we can estimate 
which value of a should lead to an amount of diffusion which is 
consistent with the above discussion. Assuming that the dispersion 
in asteroid properties as a function of heliocentric distance is con- 
sistent with diffusion generating cr(Aa/a) = O.I6(= 0.5/3) o ver a 
putative solar nebula life-time of 5 Myr ( Haisch et alJuOOlh . and 
noting that A'ao = 1.7 x 10"^^ (see Sect. l3.7.Tt . we find an upper 
limit for a ^ 5 x 10"^. Interestingly, this is very similar to the 
midplane Reynolds str ess obtained in shearing box simulations of 
discs with dead zones dFleming & Stonell2003l : iTumer et alj|2007l ; 
lllgner & NelsonluOOSl) . It would seem that our rough estimate for 
the amount of difiiision which can occur as a function of a provides 
circumstantial evidence that the solar nebula did indeed contain a 
substantial dead zone in the vicinity of the asteroid belt and giant 
planet formation region. We note, however, that caution must be 
used when interpreting this result. Confirmation of the variation of 
diffusion rate as a function of a is required using simulations which 
include vertical stratification and dead zones, and account needs to 
be made of the fact that the nebula mass changes significantly over 
a 5 Myr period. Models for the formation of giant planets require 
nebula masses approximately 3-5 times larger than the MMSN 
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Figure 24. Dispersion in the radial displacement t^x for a swarm of particles 
subject to tlie gas gravity alone. The time evolution follows a random walk 
as reasonably well described by the given relation. 



jPollack et al.l 19960 . which will increase the value of cr(Aa/a) by a 
similar factor. The disc mass will deplete over the nebula life time, 
so that the net effect of these competing factors needs to be exam- 
ined in more detail. 

Taken together, the constraints on a lead to some interesting 
conclusions regarding the strength of turbulence in the solar neb- 
ula. Shearing box simulations with dead zones generate Reynolds 
stresses at the midplane with a ~ 10"'. If this stress is accompanied 
by density waves and stochastic torques with amplitudes which fol- 
low the scaling described in Sect. |3.7.Tl then it becomes very diffi- 
cult to see how runaway growth could have caused rapid growth of 
planetary embryos via planetesimal accretion. Indeed, our results 
suggest that with a = 10"*, catastrophic disruption of 10 km sized 
bodies remains possible (although this conclusion must await firm 
confirmation from more realistic, vertically stratified simulations). 
Taken at face value, this indicates that a mechanism is required 
for forming large bodies which avoids the need to go through the 
runaw ay growth phase, such as that suggested by lohansen et al. 
( 120071) . The existence of gradients in the chemical composition of 
solar system bodies as a function of heliocentric distance, however, 
is broadly consistent with a nebula model in which there exists a 
dead zone with Reynolds stress a ^ 10"'. The existence of these 
chemical gradients however, is apparently inconsistent with a fully 
turbulent disc model with effective a » 10"'. 



3.7.3 Evolution of the eccentricity and Aa/a 

For small values of the eccentricity, the relation between the rms 
eccentricity, cr(e), and the radial velocity dispersion, cr(i',-), may be 
approximated as a-{e) ^ (t(vv)/vk, where vk is the Keplerian veloc- 
ity. In Sect. 13.6.21 we showed that for the larger planetesimals with 
ifp = 10 km, the radial velocity dispersion may be approximated as 
a random walk with cr(Vr) = Ca-(v,-) Vt, where a-(v,-) is measured in 
units of the local sound speed, Cj, and the Ccr(vr) values are listed in 
table[T] Similarly, we have shown in Sect. l3.7.Tl that the rms radial 
diffusion cr(Aa/a) = Ca-(Aa) y/t, where the Co-(Aa) values are also 
listed in table [T] Thus we expect the ratio of eccentricity changes 
to semimajor axis changes, cr(e)/cr(Aa/a) = (H/r)Ca-(v,-)/Ca-{Aa). 
Taking the value H/r = 0.05 for the models Gl - G4, we find that 
0.34 < cr{e)/cr(Aa/a) < 0.57, so that changes in the relative semi- 
major axes are similar to eccentricity changes for larger planetes- 



imals that are subject to stochastic gravitational interaction with a 
turbulent disc. 



3. 7.4 Diffusion versus type I migration of low mass protoplanets 

In addition to considering the radial drift of planetesimals due to 
gas drag, and the role of stochastic torques in potentially inhibiting 
this inward drift, we can also consider the effect that these stochas- 
tic torques might have on the type I migrati on of low mass planets 
JNelson & Papaloizod l2004l : lNelsonll2005l) . We have ascertained 
above that bodies undergoing diffusion due to stochastic torques 
may change their semimajor axes by 50 % over a disc life time 
of 5 M yr. Using the type I migration formula from iTanaka et al.l 
( 120020 for a disc model with properties which are the same as 
model Gl, we note that planets with mpi/M* = 3.5 x 10"' will 
migrate a distance equal to 50 % of their initial semimajor axis (as- 
sumed to be 5 au) via type I migration over 5 Myr. This indicates 
that direct competition between type I migration and stochastic mi- 
gration will only prevent the large scale migration of planetary bod- 
ies with masses similar to that of Mars over such nebula life times. 
It should be noted, however, that the possible radial structuring of 
protoplanetary discs by turbulence due to spatial variations in a 
(possibly due to a dead zone), may provide a means of prevent- 
ing inward drift due to th e operation of corotation torques dMassetl 
l200lUMassetetal .1120061) . 



3.8 Migration of planetesimals - local model 

The local shearing box approximation that we have adopted in this 
work does not include the effect of a radial pressure gradient in 
modifying the orbital angular velocity of the gas. Consequently, 
the gas and embedded bodies orbit at the same mean angular ve- 
locity, and planetesimals do not undergo radial migration due to 
gas drag interaction with a sub-Keplerian disc. Large planetesi- 
mals, whose interaction with the gas disc is gravity dominated, ex- 
perience stochastic gravitational forces that cause diffusion of the 
semimajor axes (or equivalently, the guiding centres of the parti- 
cle epicycles). Smaller bodies, whose evolution is dominated by 
gas drag, similarly diffuse radially due to the stochastic gas drag 
forces. We discuss both of these regimes below. 



3.8.1 Diffusion of gravity dominated bodies 

As with the larger planetesimals discussed in Sect. 13.71 for the 
global models, the larger planetesimals in the shearing box sim- 
ulations undergo radial diffusion due to the stochastic gravitational 
forces they experience. The deviation in the positions of the guiding 
centres of particles from their initial values is denoted as Ax (mea- 
sured in units of the local scale height H), and we plot the standard 
deviation of this quantity for the gravity-only particles (set "G") in 
Fig. |24l As with the global simulations, the evolution of Ax dis- 
plays random walk behaviour. As is illustrated by the fitted curve 
and error bounds in Fig. |24l the random walk can reasonably well 
approximated by the function 



(t(Ax) k Ca-iAx) Vf- fo 
with a single fitting constant C^ = 7.67 ± 1.28 x 10"^ 



(23) 
representing 



the strength of the stirring mechanism. 

The amplitude of this fitting constant for the G+D particle set 
is plotted using the light grey dash-dotted curve in Fig. 1251 (which 
also displays random-walk fitting parameters for the growth of the 
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Figure 25. Excitation amplitudes for the random walk behaviour shown in 
Fig. 1241 lines con'espond to Co- for the displacement (dash-dotted), the rms 
velocity of the guiding centres (solid), and the eccentricity (dashed lines). 
Units are in [//], [Cj], and [H/R] times [(In/Cl])''^-], respectively. 



radial velocity dispersion and eccentricity. These values are also 
listed in tableO. It is clear that for particles with Rp > 100 m, this 
coefficient has an almost constant value, indicating that the radial 
diffusion of these larger sized planetesimals is dominated by gravi- 
tational interaction with the disc. The transition region is rather nar- 
row, and G+D particles quickly approach the values corresponding 
to the gravity only particle set. 

We now compare the results shown in Fig. [24] with those ob- 
tained for the global simulations. The deviation of the position of 
the particle guiding centres is measured in units of the local scale 
height, H. We see from Fig. [24] that the best fit reaches a value 
Ax = 0.107/f after a time / = 200 orbits, and would reach a value 
Ax = 0.262/f after 1200 orbits. Dividing by the initial radial lo- 
cation of the particles, we obtain Ax/r = 0.262H/r. The shearing 
box runs have a value H/r = 0.075, when we consider that they 
are located at an orbital radius r = 5 au, such that Ax/r = 0.02. 
This compares reasonably well with the results of the global simu- 
lations, where we obtain a value of Aa/a = 0.03 for model Gl at 
t = 1200 orbits. In terms of physical parameters, model G5 is the 
most similar to the shearing box simulation, but as may be seen in 
Fig-EB the rate of radial diffusion for models Gl and G5 are very 
similar. 

As discussed in Sect. 13.71 the diffusion coefficient associated 
with the random walk in radius of the particles (or equivalently the 
diffusion of angular momentum) Dj = cr^Tcorr, where (Tt is the 
standard deviation of the fluctuating torque and Tco,r is the torque 
correlation time, averaged over the ensemble. The diffusion coeffi- 
cient found in the global simulations was found to be in reasonable 
agreement with the degree of particle diffusion obtained in those 
models. For the shearing box model, the value of ctj- is given in 
Fig. [7] and is found to be almost identical to that obtained in the 
global model Gl. The torque correlation time is given by Fig.|9] 



and has a value t,.. 



0.32 orbits. We thus see that the value of 



Ax/r = 0.02 obtained by extrapolating Fig. [24]to a time t = 1200 
orbits is in very good agreement with what would be expected 



3.8.2 Diffusion of gas drag dominated bodies 

As with the gravity dominated particles, the gas drag dominated 
particles diffuse radially over time. The radial displacement. Ax 
evolves according to a random walk, and can be reasonably well fit 
by Ea.(l23t. The fitting coefficients, Ca-iAx), are plotted in Fig. 1251 
as a function of planetesimal size. In table[2]we furthermore com- 
pile a collection of values for reference. We observe that the dif- 
fusion via drag forces roughly scales with stopping time like t^' 
for T, <c 1. When the stopping time approaches the dynamical 
time Q.', which is the case for boulders of a few metres in ra- 
dius, the dependence on t, becomes weaker because the particles 
are tightly coupled to the turbulent motion of the flow, and essen- 
tially behave like massless tracers. In this limit we see that, and 
Co-(Ax) approaches the value of the tracer particles in Fig. [25] 

To interpret the above findings in a more general context and 
be able to scale the obtained values with respect to model param- 
eters, we now consider the Schmidt number, denoted Sc = v/D|j 
This is the ratio of the momentum diffusion rate and the mass dif- 
fusion rate, and in our case the momentum diffusion v is given by 
the turbulent Maxwell and Reynolds stresses as quantified by the 
dimensionless a parameter. For the mass diffusion, we distinguish 
between Lagrangian fluid elements (represented by massless tracer 
particles) for which we obtain a diffusion coefficient Dg, and plan- 
etesimals with inertia, to which we assign the notation Dp. 

Lagrangian tracers: Before looking at the diffusion of the plan- 
etesimals themselves, we need to provide a reference value. This is 
given by SCg = v/Du, i.e., the Schmidt number for Lagrangian fluid 
elements. To determine the diffusion coeffi cient Dg, we closely fol- 



low th e approach described in Sect. 4 of jFromang & Papaloizoul 
( l2006r) . This means, we compute the velocity autocorrelation func- 
tion 



5'^(t) = <v,(z(zo,t),t)v,(zo,0)> 



(24) 



where the dependence of v, on z(zo, t) implies that we are consider- 
ing the correlation with respect to a Lagrangian fluid element. Note 
that JFromang & Papaloizoul have approximated this by v,(zo, t), 
i.e., the Eulerian velocity components (cf. their Sect. 2.3). Because 
we include the evolution of massless tracer particles, we do not rely 
on this approximation and can directly compute EQ.ll24t along the 
particle trajectories. 

The diffusion of a fluid element along the radial direction can 
then be obtained by time-integrating the corresponding component 
of the correlation tensor: 



DJr) 



f 

Jo 



5"(T')dT'. 



(25) 



from the diffusion coefficient Dj 
Ax/r = 0.021. 



cr^Tcon, the prediction being 



As discussed in lFromang & Papaloizoul (|200g), for large t this can 
be approximated by the product of the square of the rms velocity 
(vl) = S^'^XO) and the correlation time t^. of the turbulence. This 
estimate is based on the assumption that the correlations in velocity 
decay like e"'"''"'^ . To check whether this assumption is warranted, 
we plot 5g(T) and the derived D^(t). This is done separately for the 
X and y directions in the left and right panel of Fig. 1261 

Looking at the two "hockey stick" shaped curves, we see that 
the e"'"''"^ law (dashed lines) is met well for the azimuthal com- 
ponent (right panel), whereas there is considerable tension for the 



* To avoid further confusion, we adh ere to the conventions introduced in 
Sect. 2.1 of JYoudin & Lithwickl JJOOtI) . 
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1km 



CAvf) [10-3 c,] 
C^e) [10-3 ff/R] 



0.77±0.13 11.4±0.02 0.95±0.11 0.94±0.10 7.71±0.03 0.44±0.12 
3.53±0.03 45.9±3.75 2.54±0.36 2.79±0.04 45.1±3.10 0.46±0.01 
2.77±0.53 — 2.57±0.04 2.58±0.01 — 0.85±0.01 



0.01±0.18 56.5±0.01 
0.03±0.00 — 

0.03±0.00 — 



Table 2. Random walk-amplitudes Co- describing the stochastic excitation (as depicted in Fig. 1241 of the dispersion within an ensemble. Values ai'e obtained 
for the displacement Ajc, the rms radial velocity dispersion, v^, and the eccentricity e. With the exception of the first column, the applicability of the functional 
form Ccr(t - fo)''^ is of course limited to the time interval before saturation is reached (see Figs. ll3l and ll4t . 
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Figure 26. Velocity ACFs Sf (left panel), and 
massless tracers. From the ACFs {black lines, Ihs axes), we obtain time- 
integrated diffusion coefficients Dg{t) {grey line, rhs axes). The fits Sg ~ 
e-'^''^' {dashed black) result in corresponding diffusion profiles {dashed 
grey) which agree for rSlrr/il. For r 2> 1 , Dg can be approximated by 
the given theoretical estimate {dash-dotted line). 



radial direction (left panel). Note that the time domain has loga- 
rithmic scaling, which implies that the deviation in the shape of 
the curve cannot be accounted for by changing the fit parameter t^. 
Tentatively, an exponential decay law can be restored via transfor- 
mation to a generalised time variable f — » t" with a == 0.7 - we 
however leave it open whether the observed issue bears any rele- 
vance at the current stage of modelling. 

Our fitted correlation times are t*''' = 0.92 Q.' = 0.15x2;rn-' 
and r!;'' = 1.12n-' = 0.18 x 2;frn-^ The obtained numbers 
are in fact very close to the value of 0.15 orbits, reported in 
iFromang & Papaloizo ul 1120061). This supports the notion that t^ is 
characteristic for MRI-typical Mach numbers and only weakly de- 
pends on the actual amplitude of the turbulent stresses. 

The relation from Eq. ( 1251 ) is illustrated by solid grey lines in 
Fig-HD which represent the diffusion "constant" D^ as a function of 
time. For comparison, we also plot the diffusion profile correspond- 
ing to the fitted solution (dashed grey lines). These curves reason- 
ably approximate the t dependence of D„ for t :S 2n/n. Note, how- 
ever, that the fitted amplitude is somewhat smaller than the rms 
velocities and the curves saturate at a lower value. This means that 
the real correlation functions have a (stochastic) tail that leads to a 
further growth in Dg for r S 2n/Q.. 

For times sufficiently greater than the correlation time of the 
flow, the diffusion coefficient approaches a constant value approx- 
imated by Tc ( v^ ) ^ 0.035 as indicated by dashed-dotted lines in 
the two panels of Fig. 1261 In acc ordance with a recent study by 
iMadarassv & Brandenburd ( 120091) . the observed anisotropy in the 
stream-wise and cross-stream directions is rather weak. 

Translated into a Schmidt number, we yield a value of 
SCg ^ 1.6, which is considerably lower than the ratio between 
the Maxwell and Reynolds stress, which fluctuates between 3 
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Figure 27. Velocity autocorrelation function 5g(r) (solid black line) for 
particles with radius i? = 5 m. The exponential decay is now modulated by 
epicyclic oscillations and can be fitted via a function ~ cos(2;rTn-' ) e-'"''"'^ 
(white dashed line). The envelope {dotted lines) is assumed to be relevant 
for the diffusion of the guiding centre. 



and 4. Our value is a factor of about two smaller than the Eule- 
rian v alue for vertical diffusion found by iFromang & Papaloizoul 
(2006). Schmidt numbers o btained for small particles are usually 
inferr ed to be of order unity djohansen & Klaluil2005l : lTumer et al.l 
l2006r) . and our results are compatible with this. 

A Schmidt number of order unity implies that small dust 
grains that are strongly coupled to the gas will diffuse over large 
distances during protostellar disc life times. Observations of crys- 
talline silicates embedded in circumstellar discs at significant dis- 
tances from their host stars, where temperatures are too low to ex- 
plain in-situ crystallisation of amorphous grains dvan Boekel et al.l 
2004), suggest that turbulent diffusion may indeed be responsible 
for transporting grains from the hot inner regions of discs to the 
cooler outer regions. It remains to be demonstrated, however, in a 
global, turbulent disc model, whether or not a substantial number 
of grains can be transported outward against the net inward mass 
accretion flow onto the star. 

Particles with inertia: Unlike massless tracers, real particles are 
subject to the additional body forces in the rotating coordinate 
frame. Consequently, this leads to the excitation of epicyclic oscil- 
lations. In the case of Keplerian rotation, the associated angular fre- 
quency K is identical with the local rotation rate Q.. The autocorre- 
lation function for a perfect epicyclic oscillation is then simply pro- 
portional to cos(2:rrTn-'). This means that velocities are maximally 
correlated if they are apart by a full period, and anti-correlated in 
between. 

In general, the particles' motion will however deviate from 
a perfect epicycle - or, in other words, the defining elements of 
the motion will change stochastically. With both the eccentricity 
and the position of the guiding centre fluctuating, the coherence 
of the motion is lost for larger and larger times. This is illustrated 
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Table 3. Time scales and Schmidt numbers Scp = ffss c^H/Dp for the dif- 
fusion of the guiding centres of heavy particles on excited epicycHc orbits. 
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Figure 28. Attempt to determine the Schmidt number Sc = D„/Dp for the 
diffusion of the guiding centres. The diffusion coefficient is estimated via 
Df - Tc{ Vg;. ), where Vgc is obtained by smoothing over the epicycles {solid 
lines). An alternative estimate {open squares) is derived independently via 
the coefHcients C^(A,v) in Tab. [2]- note that these do not depend on any 
form of smoothing. 



in Fig. [27] where we plot the velocity ACF for particles of size 
Rp = 5 m. The epicyclic nature of the particle motion is clearly 
seen in the sinusoidal modulation of 5p. If integrated over time, the 
"diffusion" process would rather resemble a shaking motion than a 
random walk. 

As a result, this renders the direct integration approach illus- 
trated in Fig. [26] difficult. Taking the excellent agreement with the 
functional dependence ~ cos(2;rTn"')e"'"''"S we however conjec- 
ture that the diffusive part of the overall motion can be recovered 
from the exponential envelope function (dotted line in Fig.|27|(. In 
Tab. [S] we report this correlation time as a function of particle ra- 
dius Rp and comparqj it to the stopping time t^. The ratio t^/Ts is 
about unity for Rp = 5m and roughly scales with the square-root of 
the radius. 

Given we can obtain a good estimate for the rms velocity of 
the guiding centres, we can then compute Dp ^ t^ i'^'lc)- Look- 
ing at Fig.[TT] it seems less than straightforward to deconvolve the 
random motion of the guiding centre from the overall motion, espe- 
cially since the amplitude of the epicycles changes on comparable 
time scales. As a first attempt, we simply smooth-out the oscilla- 
tions applying a box-car filter in time. The thus derived Schmidt 
number Sc = D^/Dp as a function of particle size is plotted in 
Fig-HUfor various values of the filter size. 

If we omit the smoothing, the results (although not strictly 
correct) reflect the kinetic energy in the epicyclic mode itself. In 
this case, the gradual decoupling of the particles implies that { v^ > 
decreases while t^ increases accordingly. Interestingly, this leaves 
the overall particle diffusion rate unchanged. If we increase the 



^ The ratio Ts /tc should not be confused with the Stokes number St = Ts/t^ 
which measures the particle coupling with respect to the eddy velocity. 



filter size, the epicyclic motion is attenuated and the size depen- 
dence gradually becomes steeper. Reasonable convergence is ob- 
tained when the filter size approaches the critical value, equal to 
the epicyclic period. In Tab.[3]we list the corresponding numbers. 

Note that already for metre-sized objects the particle diffusiv- 
ity is reduced relative to Lagrangian tracers. We attribute this to the 
effect that part of the kinetic energy is absorbed in ordered epicyclic 
m otion and is thus not availa ble for particle diffusion (cf. Sect. 3.2 
in lYoudin & Lithwickll2007h . The partial slippage of the particles 
through the gas as a result of weaker drag acceleration also causes 
the particle guiding centre velocities to be reduced relative to the 
fluid turbulent velocities. 

It might be seen as a deficiency of the approach that it depends 
on the proper choice of a box-car filter function. To reinforce the 
validity of the derived Schmidt numbers, we therefore provide an 
independent estimate from the squares of the coefficients Co-(Ajc) 
from Fig. [25] Defining the growing spread in the particles position 
with time, these can similarly be interpreted as a diffusion constant. 
For consistency, we have also checked that the fit constants Ca-(Ax) 
do not change significantly, when the epicyclic part of the motion 
is smoothed. Considering that this second estimate solely relies on 
positions while the original one is derived from velocity correla- 
tions, the agreement seen in Fig. [28] (where we plot the different 
estimates for the Schmidt number) is quite striking. 



4 CONCLUSIONS 

We have presented results from both local and global MHD simu- 
lations which examine the evolution of planetesimals and boulders 
of different size embedded in turbulent protoplanetary disc models. 
The main aims of this work are to: determine the magnitude of the 
velocity dispersion which is excited in a planetesimal swarm; de- 
termine the rate of radial diffusion of the planetesimals; determine 
under which conditions broad agreement is obtained in the results 
from local shearing box simulations and global disc simulations. 

We find that the magnitude of density fluctuations, and the as- 
sociated stochastic forces experienced by embedded planetesimals 
is sensitive to the dimensions of shearing box models. This appears 
to be related to the fact that accurate modelling of the excitation and 
non-linear steepening of spiral density waves requires boxes which 
are elo ngated in the azimuthal direction i JHeinemann & Papaloizoui 
l2009a) . The two-point correlation function for the perturbed sur- 
face density also shows that coherent structures stretched by a dis- 
tance ^ 6H are a common feature of MHD turbulence in discs, 
indicating that shearing boxes in excess of this length scale are re- 
quired to prevent the premature truncation of the gravitational field 
which results from these structures. In addition to the magnitude 
of the stochastic torques being dependent on the box size, we also 
found that the correlation time of the fluctuating torques also de- 
pends on the box size and aspect ratio, with smaller boxes gener- 
ating shorter correlation times. This is apparently due to the ability 
of waves to propagate across the box on multiple occasions due 
to the periodic boundary conditions employed, combined with the 
rate at which they shear past the planet due to the background flow. 
Shearing boxes with dimensions 4H x 16H x 2H were found to 
provide converged results which agree well with the results from 
global simulations. 

We find that both global simulations and local shearing box 
simulations predict that rapid excitation of planetesimal random 
velocities is expected in fully turbulent disc models whose local 
surface densities are similar to the minimum mass solar nebula. A 



© 2002 RAS, MNRAS 000.[Tll24l 



22 Nelson & Gressel 



model whose turbulent stresses generate a = 0.035 leads to rapid 
growth of the radial velocity dispersion, cr(Vr), via a random walk, 
such that after 1200 orbits at 5 au (r(vr) = 200 m s"' . A model with 
weaker turbulence (a = 0.017) gave rise to (t(v,) = 166ms"' 
over the same evolution time. These random velocities are much 
larger than either the escape velocities from planetesimals with 
sizes 1km or 10 km, or the catastrophic disruption thresholds for 
collisions between b odies of similar size JBenz & Asphauall999l : 
IStewart & Leinhardt 2009), suggesting that planetesimal collisions 
occurring in fully turbulent discs similar to those considered in this 
paper will result in destruction of the planetesimals. We find that the 
expected equilibrium velocity dispersion for 10 km sized planetes- 
imals scales weakly with the turbulent stress parameter such that 
cr(Vr) = Ky^a"-^. Extrapolation of the simulation results presented 
in this paper to low values of a suggest that even with a = 10"*, 
bodies with size 10 km are likely to have random velocities which 
will cause catastrophic disruption during mutual collisions. Values 
of cr(v,.) small enough to allow runaway growth of planetesimals 
to occur will be even more difficult to achieve. For smaller bodies, 
which are more tightly coupled to the gas via gas drag, we find that 
the equilibrium velocity dispersion decreases, and reaches a min- 
imum of cr(Vr) - 20 ms"' for bodies of size 50 m in a disc with 
a = 0.035. For smaller bodies than these the increasing influence 
of gas drag causes the velocity dispersion to increase with decreas- 
ing size, and boulders with Rp = Im develop random velocities 
which closely match the turbulent gas motions. 

In addition to driving the growth of the velocity dispersion, or 
equivalently the eccentricity, the stochastic forces experienced by 
the planetesimals cause them to diffuse radially through the disc. 
For large planetesimals with Rp = 10 km, radial drift through the 
disc due to gas drag occurs very slowly, and their radial motion 
is expected to be dominated by diffusion caused by the stochas- 
tic gravitational field of the turbulent disc. Indeed, after a putative 
disc life time of 5 Myr, the rms relative change in semimajor axis, 
cr(Aa/a) for a planetesimal swarm located at 5 au in a disc with 
a = 0.035 will be 50%, which would have been sufficient to cause 
large scale migration of the small body populations of the Solar 
System during its formation. The fact that such large scale migra- 
tion of large planetesimals does not appear to have happened allows 
us to put constraints on the strength of a in the solar nebula. A value 
of Q- ^ 10"^ would lead to a more reasonable cr(Aa/a) ==0.1, con- 
sistent with the notion that the solar nebula had a dead zone in the 
region where planet formation took place. 

Smaller planetesimals with Rp = 10 - 100 m are expected to 
undergo large-scale radial drift due to gas drag, and this is observed 
in our simulations. The migration rate of 10 m size bodies in the 
turbulent discs, relative to that expected for an equivalent laminar 
disc model, is found to be very similar. Migration in turbulent discs 
is found to be slower by a factor of approximately 2-3, and this 
is due to radial structuring of the disc's density profile, creating 
regions where the pressure gradient (and hence rotational velocity 
profile) is modified. 

There have been a number of previous studies of the dynam- 
ical^ evolution of planetary bodies embedded in turbulent discs. 
[Nelson & Papaloizou ( 2004) examined the torques experienced by 
low mass planets and suggested that random walk behaviour for 
suc h bodies should be expected, as has been observed in this pa- 
per. iNelsonI ( 120051) examined the orbital evolution of planets em- 
bedded in turbulent discs, with particular emphasis on the semima- 
jor axis and eccentricity changes. The strength of the turbulence 
in that model was somewhat weaker than considered in this pa- 
per (a - 5 X 10"^), although the disc model was somewhat more 



massive, and in general it was found that the response of embed- 
ded bodies to the stochastic forcing was stronger in that study, 
with larger changes in semimajor axis and eccentricity being ob- 
served. The primary reason for this discrepancy appears to be the 
fac t that a persistent vortex formed in the disc model considered 
bv lNelsonI (;2005). Su c h vor tices were also reported in the work 
of iFromang & Nelsonl ( 120051) ■ and interaction between embedded 
bodies and such structures can lead to significant modification of 
the semimajor axes and excitation of eccentricity. The models pre- 
sented in this paper do not have such vortex features, but during 
the early stages of this project disc models with lower values of a 
were found to generate such flow features quite readily. One pos- 
sibility is that global disc models which are close to the limits of 
resolving the MRI generate substantial variations in o- as a function 
of radius, which in turn generate localised pressure maxima which 
are particularly pron e to the formation of vortices cHawlevll 19871 : 
iLovelace et alj|l999f) . Models with smaller values of the plasma /3 
parameter have stronger field strengths and thus resolve the MRI 
more easily, and in addition the stronger fields may act to inhibit 
vortex formation through the action of local magnetic stresses. It is 
thes e models that we h ave considered in this paper. 

lYang et alj ( 120091) recently considered the evolution of plan- 
etesimal swarms using non stratified shearing box simulations. The 
planetesimals in that study experienced the fluctuating gravitational 
field of the disc, but did not experience gas drag forces. Although 
random walk behaviour of the particles in that study was observed, 
in general it was found to be significantly weaker than we report 
here, and one conclusion reached by these authors is that turbu- 
lence is unlikely to cause the catastrophic disruption of planetesi- 
mals. The primary reason for this appears to be the fact that most 
simulations performed by lYang et al.l ( 120091) used shearing boxes 
with 2H X 2H x 2H. Test calculations with larger boxes presented 
in that paper indicated a strong increase in the response of the par- 
ticles to the turbulence forcing as a function of increasing box size, 
suggesting that the main reason for the discrepancy in our results is 
due to this effect. 

Ilda et al.l j2008l) recently considered the evolution of planetes- 
imals embedded in turbulent protoplanetary discs by means of N- 
body simulations combined with a prescription for planetesimal 
stirring based on the work of lLaughlin et alj ( 12004) . In basic agree- 
ment with the results we have presented in this paper, they show 
that turbulence leads to the excitation of a large velocity disper- 
sion, which is likely to cause catastrophic disruption of planetes- 
imals rather than growth following mutual collisions, for a wide 
range of turbulent strengths. 

The simulations we have presented here use the simplest pos- 
sible numerical set-up: ideal MHD in non stratified disc models. In 
a future paper we will present the results of a similar study using 
vertically stratified models with and without dead zones. This fu- 
ture study will tell us whether we need to examine new paradigms 
for the rapid growth of planetesimals, or whether instead runaway 
growth can indeed occur within a dead zone of an otherwise turbu- 
lent protoplanetary disc. 
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